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1.  INTRODUCTION 


For  numerous  years,  personnel  of  the  Battlefield  Environment  (BE)  Directorate  of  the  U.S. 
Army  Research  Laboratory  (formerly  the  U.S.  Army  Atmospheric  Sciences  Laboratory 
(ASL)  of  Laboratory  Command)  have  been  actively  investigating  atmospheric  influences  on 
Army  target  acquisition  systems.  Among  the  research  topics  impacting  target  acquisition 
are  the  areas  of  atmospheric  image  distortion  due  to  turbulence,  optical  refraction,  acoustic 
propagation,  scattering  of  acoustic  waves  by  turbulence,  and  turbulence  effects  on  laser 
propagation.  These  topics  have  been  considered  at  length  and  continue  to  receive  attention. 

The  nature  of  this  group  of  problems  is  such  that  the  atmosphere's  vertical  structure 
near  the  surface  plays  a  critical  role  in  the  performance  of  Army  systems.  With  optical 
refraction,  the  vertical  structure  of  temperature  is  important.  For  electromagnetic  energy 
propagation,  either  imaging  or  laser  propagation,  the  critical  parameters  are  the  inner  scale 
of  turbulence,  the  outer  scale  of  turbulence,  L^,  and  the  index  of  refraction  structure 
parameter,  C^.  is  composed  of  a  temperature  fluctuation  structure  parameter  (C^),  a 
related  humidity  term,  a  cross  correlation  term  characterizing  interactions  between  tern* 
peratur<‘  an<l  humidity  fluctuations,  and  a  pressure  term.*  For  acoustic  wave  propagation 
the  tcmpc'rature,  pr<\sstire,  and  wind  vertical  structures  are  important.  For  scattering 
of  acoustic  waves,  the  temperature  (C^)  and  wind  (Cy)  structure  parameters  must  be 
considered. 

1.1  History 

Since  these  propagation  problems  all  have  the  vertical  structiure  of  the  near  surface 
atmosphere  in  common,  it  would  be  useful  to  be  able  to  predict  the  vertical  structure 
of  this  region.  Recently,  H.  Rachele  and  A.  Tunick  of  BE  have  taken  up  this  noble 
goal.  Several  papers  have  been  produced  concerning  new  techniques  for  treating  sensible 
and  latent  heat  flux  effects  on  C^.  The  applications  of  the  Rachele/ Tunick  work  have 
been  in  the  eurea  of  turbulence  effects  on  image  propagation.  Similarly,  in  the  early  80’s, 
K.  Kunkel  and  D.  Walters  of  ASL  engaged  in  modeling  the  near  surface  atmosphere  for 
charactcrizatiiJii  of  atm<)S])h<;ric  influences  on  high-energy  laser  propagation. 

The  work  presented  in  this  n'port  wjus  accompli.shcd  between  1982  and  1987  while  the 
author  wiis  working  on  chiU'Hct<'ri/.Htion  of  the  atmosphere  to  predict  the  probability 
density  function  of  atmospheric  refractive  effects  in  desert  locations.  The  work  is  therefore 
somewhat  dated,  although  the  treatment  of  foliated  layers  is  believed  to  be  unique.  It  is 
hoped  this  report  will  provide  a  fixed  point  of  reference  concerning  the  approach  taken  at 
the  time  to  handle  the  surface  energy  budget  problem.  The  works  of  Rachele  and  Tunick 
did  not  exist  at  the  time  of  the  work,  and  no  attempt  has  been  made  to  incorporate  their 
findings  in  this  report.  As  such,  for  a  wet  enviromnent  the  model  presented  will  be  less 
than  optimal,  yet  for  its  intended  purpose  (refraction  in  deserts)  the  model  should  be 
appropriate. 


*Usually  the  pressure  term  is  ignored,  as  will  be  done  here. 
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1.2  Background 


The  goal  of  a  surface  energy  budget  model  is  to  simulate  the  heat  fluxes  at  the  earth’s 
surface  and,  in  so  doing,  to  recreate  the  environment  wherein  the  vertical  structures  of  the 
atmosphere  exist.  This  approach  assumes  that  there  exist  speciflc  linkages  between  the 
fluxes  at  the  surface  and  the  vertical  structvires.  This  linkage  is  handled  through  similarity 
theory  and  the  flux-profile  relationships. 

Since  Halstead  et  al.  (1957)  developed  an  analog  computer  to  simulate  the  surface 
energy  budget,  many  authors  have  attempted  to  derive  equations  that  accurately  calculate 
turbulent  heat  fluxes  in  the  surface  layer.  When  this  research  is  coupled  to  flux-proflle 
theory  (Dyer,  1974),  atmospheric  structvire  proflle  predictions  for  near-surface  conditions 
are  possible. 


2.  MODEL  DISTINCTIVES  AND  THEORY 

In  general,  the  theory  discussed  here  deals  with  the  basic  equations  necessary  to  calculate 
the  surface  layer  aerodynamic  parameters  of  friction  velocity,  ti,,  Obukhov  length,  L 
(Obukhov,  1946),  soil  temperature,  and  foliage  temperature.  Once  these  parameters  are 
known,  only  a  single  step  is  needed  to  produce  complete  estimates  of  the  vertical  structures 
mentioned  in  secticwi  1.  To  compute  u,  and  L,  the  sensible  heat  flux,  Hsh,  must  first  be 
found.  To  find  th<!  .sensibU;  Iwat  flux,  the  force-restore  surface  energy  budget  method  is 
used.  This  method  is  simultaneously  coupled  with  an  iterative  Obukhov  length  calculation. 


2.1  Surface  Energy  Budget 

The  surface  energy  budget  approach  developed  here  is  based  on  Deardorff’s  (1978)  model 
that  included  vegetation  layer  effects.  This  model  had  been  progranuned  in  FORTRAN  by 
K.  Kunkel  of  ASL  prior  to  the  time  the  author  began  work  on  the  model.  The  Deardorff 
model  calculates  the  radiative,  convective,  and  conductive  energy  fluxes  firom  a  grotmd 
sinrface  layer  and  a  foliage  layer.  Here  the  use  of  the  term  foliage  layer  refers  to  an  equivalent 
surface  of  plants  that  obstructs  the  transit  of  radiation  flowing  from  the  atmosphere  to 
the  earth.  The  mean  height  of  this  layer  is  variable,  as  well  as  the  overall  fraction  of  the 
ground  surface  covered.  The  fractional  foliage  coverage  is  denoted  by  the  term  a/  and 
repres<uits  that  fraction  of  th«;  surface  that  intercepts/blocks  direct  sunlight  from  reaching 
the  ground.  In  this  model  the  foliage  layer  can  be  thought  of  as  a  low  lying  layer  of  weeds, 
grass,  or  bushes. 

The  fluxes  to  and  from  the  surface  may  be  divided  between  those  above  the  foliated 
layer  (denoted  by  h),  fluxes  between  the  foliage  and  the  smface  (denoted  by  g),  and  the 
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Figure  1.  Modeled  fluxes  at  soil  surface  and  modeled  foliage  layer  surface. 

conductive  heat  flux  into  the  soil  (denoted  by  G).  Figure  1  shows  the  directionality  of  the 
fluxes  involved. 

The  fluxes  are  further  divided  by  the  type  of  energy  being  transported.  There  are 
shortwave  (5)  and  longwave  {Rl)  radiative  fluxes,*  convective  fluxes  in  latent  form  (due 
to  evaporation  or  condensation)  {L„E)  and  in  “sensible”  form  (due  to  warming  or  cooling 
of  air  at  the  surface)  (/f.v),  and  the  conductive  flux  (G)  of  heat  into  the  soil.  Radiative 
fluxes  can  be  directed  either  upward  or  downward,  while  sensible  and  latent  heats  are  only 
handled  by  single  functions  (having  positive  senses  when  the  fluxes  are  directed  away  from 
the  surface). 

2.1.1  Controlling  Variables 

Excluding  the  shortwave  fluxes,  the  remaining  fluxes  depend  on  one  wind,  six  temperature, 
and  four  humidity  variables.  Time  variations  in  these  variables  to  a  great  degree  determine 
the  behavior  of  the  vertical  structures  to  be  predicted.  The  wind  variable,  u,  is  a  measured 


*Shortwave  usually  refers  to  all  radiation  at  wavelengths  shorter  than  4  /xm.  Longwave 
refers  to  all  radiation  at  wavelengths  longer  than  4  fita.  The  separation  occurs  naturally 
since  solar  ratliation  is  negligible  at  wavelengths  greater  than  4  pm,  while  terrestrial 
(graybody)  raxliation  has  negligible  contributions  at  wavelengths  below  4  pm.  (See  Oke 
(1978),  page  11.) 
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input  to  the  model.  It  is  measured  at  a  model  input  reference  height,  z,.  The  first  five 
of  the  temperature  variables  are  Bg,  the  temperattire  of  the  ground  at  the  stirface;  $2,  the 
temperature  of  the  groimd  a  characteristic  distance  below  the  surface,  where  this  distance 
is  dependent  on  soil  type;  Bf,  the  foliage  temperature;  B^,  the  measured  air  temperatvu% 
at  the  reference  height;  and  a  simulated  air  layer  temperature  at  the  reference  height 
that  is  driven  by  the  surface  heat  fiuxes  and  a  restoring  term  based  on  B^.* 

The  groimd  temperature  is  modeled  as  a  function  of  the  driving  energy  fiuxes  at  the  surface, 
as  modulated  by  what  Hoifert  and  Storch  (1979)  call  the  thermal  inertia  of  the  ground. 
The  thermal  inertia  depends  on  the  depth  (di )  to  which  the  thermal  wave  can  penetrate 
the  soil  in  a  day.  This  penetration  depth  is  modeled  as 


dj  —  y/H»  ,  (1) 

where  k,  is  the  soil’s  thermal  diffusivity  in  m*s“*,  and  is  the  number  of  seconds  in  a 
day  (86400  s).  The  deep  ground  temperature  62  is  modeled  as  controlled  by  the  depth  of 
penetration  of  the  thermal  wave  on  a  yearly  basis  (^2  =  v^365di  =  19.1  di). 

The  sixth  temperature  variablt  is  the  upper-air  temperature  (^,a),  which  is  used  in  the 
estimation  of  the  downward  longwave  sky  radiation  from  clouds.  The  value  of  B^a  is  based 
on  the  measured  air  temperature  and  the  type  of  clouds  present.  A  simple  estimation 
scheme  is  used;  for  low  clouds  B^a  «  —  8,  for  medium  clouds  ^*0  ^  B^  —  16,  and  for 

high  clouds  B^a  ^  Ba  —  24.  The  predictions  are  based  on  assumed  heights  of  cloud  bases 
of  about  1200,  2400,  and  3600  m,  respectively,  using  a  6.5  K/km  tempierature  lapse  rate 
(the  moist  adiabatic  lapse  rate). 

The  4  humidity  variables  used  in  the  model  are  the  atmospheric  relative  humidity  measured 
at  the  station  height  {Ru)^  the  amount  of  dew  (in  kg/m^)  on  the  foliage  (tW4eu>)»  w»d  the 
fractional  water  contents  (by  volume)  within  the  first  10  cm  of  soil  (iWj)  and  within  the 
first  50  cm  of  soil  (1^2)-  The  distances  10  cm  and  50  cm  roughly  equate  to  the  depths 
of  soil  affected  by  the  daily  {d\ )  and  yearly  (1/2)  variation  of  temperature.  Thus  the  soil 
characteristics  in  the  layer  nearest  the  surface  will  be  characterized  by  Bg  and  Wg.  Similarly, 
B2  and  u>2  will  be  associated  with  the  deeper  layer.  The  importance  of  this  association  is 
made  clear  in  section  2.3  when  the  equations  for  soil  thermal  characteristics,  which  depend 
on  the  moisture  content,  are  discussed. 

Below  the  depth  of  50  cm,  the  model  assumes  the  ground  has  no  fluctuations  in  temper¬ 
ature.  Oke  (1978)  shows  a  yearly  peak-to-peak  thermal  wave  amplitude  of  approximately 


*In  all  equations  h<!reafter  B'^  will  be  rejilaced  by  Bt,  which  is  considered  descriptive  of  the 
air  teuip<;rature  at  height  z,.  In  general  Bt  will  not  equal  the  measured  air  temperature 
since  it  is  assumed  th<*  surface  temperature  controls  the  variation  of  air  temperature  near 
the  surface  rather  than  air  temperature  driving  surface  temperature  variations.  B^  is 
the  means  of  handling  this  driving  mechanism.  A  combination  of  surface  flux  driven 
temperature  variations  and  a  stabilizing  term  allow  the  model  to  pace  the  behavior  of  the 
near  surface  atmosphere.  The  stabilizing  term,  moreover,  avoids  the  model  phenomenon  of 
surface  temperature  “creep”  that  often  plagues  long  term  prediction  models.  This  change 
was  effected  after  Kunkel’s  (1985)  review  of  the  model  and  upon  his  suggestion. 
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10  K  at  1  m  below  the  surface.  Thus  the  approximation  used  in  the  model  could 
be  improved,  though  it  is  unlikely  to  strongly  aifect  behavior  of  atmospheric  structure 
predictions  since  the  heat  exchange  with  the  soil  below  50  cm  is  small  in  comparison  to 
the  radiative  nd  convective  fluxes  at  the  surface. 


2.1.2  Introduction  to  Heat  Flux  Equations 

The  temperatures  0g  and  6/  can  be  found  as  a  function  of  time  by  solving  the  surface 
energy  budget  equations.  The  equation  for  the  surface  energy  budget  at  the  foliated  layer 
can,  in  fact,  be  cast  as  a  direct  function  of  the  foliage  temperature.  Therefore  solution 
of  one  of  the  energy  budget  equations  allows  direct  estimation  of  one  of  the  temperature 
parameters  of  interest. 

Since  the  foliated  layer’s  ability  to  store  heat  is  assumed  negligibly  small,  the  fluxes  arriving 
at  the  foliage  layer  from  above  and  below  must  sum  to  zero.  FVom  figure  1  this  sum  Ccin 
be  expressed  as 

Shi  +  Hsg  +  LEg  = 

Shi  H-Lhi  +  H-Lgi  +  Hsh  +  LEh. 

The  foliage  temperature  is  found  by  means  of  a  Newton- Raphson  method  (see  section  2.2), 
where  the  foliage  temperature  defines  a  unique  solution  to  equation  (2). 

The  condition  that  mjikes  the  calculation  of  energy  fluxes  at  the  foliated  layer  critically 
important  is  precisely  because  the  foliage  surface  has  no  thermal  inertia  and  because  the 
leaves  have  more  surface  area  than  the  ground  below.  That  is,  while  the  fraction  of  direct 
sunlight  blocked  by  the  leaves  rarely  exceeds  70  percent  (Deardorff,  1978),  the  exposed 
leaf  surface  area  (available  for  evapotranspiration  and  sensible  heat  exchtinge)  is  about 
seven  times  the  surface  area  of  terrain  overshadowed.  Thus  a  fully  vegetated  surface  can 
b(;  much  more  of  an  influence  on  the  convective  fluxes  of  momentum,  sensible  heat,  and 
latent  heat  than  the  underlying  surface.  Foliated  layer  effects  therefore  have  a  significant 
impact  on  atmospheric  structure  predictions,  which  will  alter  propagation  model  results, 
which  will  in  turn  influence  system  performance  under  those  atmospheric  conditions.  To 
properly  model  foliage  effects,  a  surface  energy  budget  model  was  modified  by  Deardorff 
to  account  for  the  interactions  at  the  foliated  layer. 

The  other  key  temperature  to  be  predicted  is  the  soil  surface  temperature.  To  predict  this 
temperature  as  a  function  of  time  requires  a  different  method  firom  that  used  to  compute 
foliage  temperature  because  of  the  thermal  inertia  of  the  groimd.  Here,  heat  is  conducted 
into  and  stored  in  the  soil  as  well  as  being  exchanged  in  fluxes  above  the  surface.  These 
storage  terms  make  the  Newton-Raphson  method  inappropriate. 

The  amount  of  energy  in  the  soil  as  a  whole  is  characterized  by  the  variable  02-  Yet  it 
is  only  the  soil’s  temperature  at  the  surface  (0g)  that  is  involved  in  equations  describing 
energy  fluxes  to  the  atmosphere  above.  Thus  two  equations  are  required,  one  describing 
tlu^  energy  fluxes  at  the  surface,  which  arc  necessarily  in  balance, 

G  =  Sgi  -  Sgi  -1-  Rig,  —  Rf.g^  —  Hsg  “  LEg,  (3) 
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and  a  second  equation  describing  the  actual  ground  surface  temperature  change, 


^9  _  £1^  .  C2(g2  -  6g) 

dt  /i  T,.,  • 

This  is  called  the  force-restore  method,  based  on  the  theories  of  Bhumralkar  (1975)  and 
Blackadar  (1976).  Blackadar  used  the  values  ci  =  3.72  and  C2  =  7.4  for  the  transfer 
coefficients,  and  /i  is  Hoffert  and  Storch’s  (1979)  thermal  inertia 


h  =  P»c,di, 


(5) 


where  p,  is  the  soil  density  and  c,  is  the  soil’s  specific  heat.  (Together,  the  quantity  p,  c, 
is  the  soil  heat  capacity.) 

Though  62  is  approximately  constant  over  the  span  of  a  few  days,  when  modeling  long 
time  spans  it  must  be  calculated  by  the  equation 


092  G 


I2  =  P»C,d2f 


d2  —  ^/365  TfinyKf  . 


(G) 


To  model  the  annual  trend  of  the  deep  soil  temperature,  R.  Johnson  (then  of  ASL) 
developed  the  following  approximation  to  be  used  to  initialize  62  for  desert  environmental 
conditions: 

62  291.66  +  9  sin  ^360-^-— 

\  366 

where  $2  is  the  temperature  given  in  degrees  Kelvin  and  Jo  is  the  Julian  date.  When  the 
Johnson  approximation  is  used,  model  behavior  of  the  G  flux  will  rapidly  approaich  realistic 
values  since  there  will  not  be  transfer  of  energy  between  the  deep  soil  amd  the  near  surface 
soil  with  the  wrong  sign  for  that  time  of  year  (when  considering  a  desert  environment  at 
approximately  30“  N  latitude).  Approximations  similar  to  the  Johnson  equation  could  be 
developed  for  other  environments  as  well. 


2.1.S  Shortwavv.  Radiation 

Each  of  the  flux  terms  in  equations  (2)  and  (3)  will  now  be  considered  in  explicit  form. 
The  first  to  be  considered  is  the  shortwave  radiation  from  the  sky,  Shi  •  This  flux  depends 
on  the  solar  zenith  angle  and  the  amoimt  of  cloud  cover.  The  amount  of  cloud  cover 
(Cc)  must  be  input  to  the  model.  This  input  is  normally  accepted  for  a  single  layer,  but 
provisions  in  the  model  allow  for  input  of  separate  skycover  fractions  at  low,  mediiun,  and 
high  cloud  heights  utilizing  the  cloud  model  of  Shapiro  (1982). 

One  of  the  key  aspects  of  any  radiation  model  is  to  locate  the  sim  within  the  sky.  Several 
additional  model  inputs  and  calculations  are  required  to  compute  the  solar  zenith  angle. 
In  subroutine  SUM  the  solar  zenith  is  computed  from  a  knowledge  of  local  standard  time, 
latitude,  longitude,  and  day  of  the  year.  The  following  procedure  is  used: 
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The  time  of  local  noon  is  given  in  minutes  aftei  midnight  as 

E  =  720-er  +  4AA,  (7) 

where  er  is  the  equation  of  time  (Smithsonian  Meteorological  Tables,  1951)  approximated 
by  (Woolf,  1968) 


-“(i 


in(A)  cos(A)  sin(2A)  cos(2A) ' 
09258  ~  233.155  6.50157  ~  16.45278. 


^  =  S’ 

where  all  arguments  to  trigonometric  functions  are  given  in  degrees. 

The  quantity  720  —  er  computes  the  time  of  solar  noon,  but  only  for  locations  at  standard 
meridians  (Aatdm)-  The  Aa  term  in  equation  (7)  is  the  means  of  correcting  for  solar  noon 
at  a  local  longitude.  Assuming  east  longitudes  are  treated  as  negative  and  west  longitudes 
as  positive,  Aa  is  defined  as 

“  Aiocal  (^1) 

whore  Alocal  is  the  local  longitude  in  degrees. 

Using  this  convention,  anytime  a  location  is  west  of  the  standard  meridian  used  to  measure 
clock  hours,  the  time  of  solar  noon  will  be  4  minutes  later  per  degree  of  difference  between 
the  local  longitude  emd  the  standeird  meridian.  Normally  this  correction  will  be  small,  but 
for  some  time  zones  it  may  be  up  to  an  hour’s  difference. 

The  second  quantity  to  be  calculated  by  the  SUH  subroutine  is  the  solar  declination,  Ds- 
A  formulation  based  on  Woolf  (1968)  is  used. 


Ds  =  sin  '(sin(23.4438)  sin(S)), 


(10a) 


S  =  A  +  279.9348  +  1.914827  sin(A)  -  (lofc) 

'  ^  12.5747  50.1555  617.284’  ^  ’ 

where  all  arguments  to  the  trigonometric  functions  are  assumed  measured  in  degrees.* 

Dt'dination  is  considtTed  constant  over  the  course  of  a  single  Julian  date. 

The  SUN  routine  then  computes  ffj,  the  fraction  of  the  day  between  sunrise  and  local  noon. 
By  calculating  this  quantity,  sunrise  and  sunset  times  can  be  determined  about  the  local 


noon. 


cos(2 


sin P  —  Sint  sin  Ds 
cost  cos  Ds 


*Appropriate  conversions  to  radians  will  be  reqmred  for  some  compilers. 


13 


where  t  is  the  local  latitude,  and  ^  (equal  to  -50  minutes  of  arc)  is  the  elevation  angle  when 
stmrise  occiu^  is  nonzero  due  to  astronomical  refraction  effects).  Sunrise  and  sunset 
times  (in  minutes  from  midnight)  are  then 


I»*nra  —  “  “  1440 


(12a) 


/«n.<  =  H-|-1440frrf. 


(125) 


Given  the  declination,  sunrise  and  simset  times,  and  the  local  time  in  minutes  past 
midnight,  one  can  compute  the  solar  energy  flux  reaching  the  top  of  the  earth’s  atmosphere 
using  a  plane  parallel  assumption.  Subroutine  CLOUD  is  called  to  make  this  calculation. 

Within  subroutine  CLOUD  (based  on  Shapiro  (1982)),  reflections  between  cloud  layers  are 
treated.  First,  the  solar  radiation  is  calculated  from  the  average  solar  radiation  intensity 
at  the  earth’s  mean  orbit  radius  (5o(m«an)=1353  W/m^)  and  is  corrected  by  factors 
accounting  for  the  actual  distance  to  the  sun.  This  value  is  then  adjusted  according 
to  the  cosine  of  the  zenith  angle,  the  effective  surface  albedo,  and  the  cloud  type  and  cloud 
amoimts  in  three  cloud  layers. 

Mathematically,  this  process  is  described  in  the  following  equations.  The  actual  solar 
irradiance  at  the  earth’s  distance  from  the  sun  is  determined  by  accounting  for  the  distance 
variations  of  the  earth  from  the  sun  as  a  function  of  day  of  the  year  through  the  equation 


5o 


^0(fn«an) 


(13a) 


where  Di.si/Dist  represents  the  ratio  of  the  average  to  the  actual  distance  to  the  sun: 


Dist 

Dist 


^1.00014  -I-  0.016726  cos 


^360(Jp-2)^^ 


(135) 


The  irradiance  at  a  plane  perpendicular  to  the  earth-sim  axis  must  now  be  translated  into 
the  energy  entering  a  plane  parallel  atmosphere  through  use  of  the  cosine  of  the  zenith 
angle.  The  cosine  of  the  zenith  angle  (cos  Q  is  determined  from  the  latitude  (£),  declination 
(Ds),  and  hour  angle  (Hr)  as 


cosC  =  sin£  sin i>5  +  cosi  cos  Ds  cosi/r, 


(14) 


where  time  t  in  minutes  past  midnight  is  converted  into  an  angle  through  Hr  =  360t/1440. 
Thus 


Si  =  So  cosC 


is  the  energy  flux  from  the  sun  (in  W/m*)  at  the  top  of  the  atmosphere. 

This  flux  is  propagated  to  the  surface  where  a  portion  of  it  is  absorbed.  The  absorbed 
portion  is  (1  —  o*),  where  is  the  effective  albedo  of  the  surface.  The  effective  albedo 
is  the  average  reflectivity  of  the  surface.  This  factor  considers  the  fractional  reflectance 
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from  that  portion  of  the  surface  covered  by  foliage  (<7/)  and  reflection  from  the  fractional 
exposed  surface  area  (1  —  a/): 

a«  =  <7/  a/  +  (1  -  <r/)*  a,  Aoo(,),  (15a) 

■Aoo(*)  =  (1  -  a,  Of  a/)“* ,  (155) 

where  a/  is  the  mean  shortwave  reflectivity  (albedo)  of  the  foliage,  and  Og  is  the  albedo  of 
the  soil.  The  effective  albedo  is  thus  composed  of  a  component  representing  reflection  from 
the  foliage  plus  a  component  representing  that  portion  of  the  incident  energy  that  passes 
through  the  foliage  layer,  reflects  off  the  surface,  and  passes  back  through  the  foliage  layer. 
The  factor  Aoo(«)  represents  the  effects  of  contributions  from  multiple  reflections  between 
the  foliage  layer  and  the  surface  before  the  energy  escapes  into  the  upward  fltix  above  the 
layer.* 

In  the  CLOUD  subroutine,  a  complicated  series  of  reflection  and  transmission  calculations 
between  different  cloud  layers  and  the  surface  is  made.  Shapiro  (1982)  gives  a  fuller 
explanation  of  this  procedure.  Here  it  will  be  expressed  simply  as 

Ski=Soxco8C,  (16) 

where  X  represents  the  fractional  attenuation  of  sunlight  due  to  the  cloud  filled  atmosphere. 

2.1.4  AtmospheTtc  “Longwave”  Radiation  Calculation 

As  previously  explained,  atmospheric  radiation  is  characterized  by  infrared  emissions  from 
atmospheric  constituents.^  In  this  model  the  downward  longwave  flux  from  the  atmosphere 
is  parameterized  using  the  equation 

+(!-€,.)  Cci  <7  (17) 

accounting  for  near-surface  molecular  emissions  and  emissions  from  clouds.  In  this  equation 
6a  is  the  average  emissivity  over  the  infrared  spectrum  of  the  air  near  the  surface  (thus 
distinguishing  it  from  blackbody  radiation  where  the  emissivity  is  imity). 

The  cloud  layer  contribution  to  the  longwave  flux  is  characterized  by  the  second  term  on 
the  right  in  equation  (17).  The  emissivity  of  the  cloud  layer  is  assumed  to  be  unity.  Cct  is 
the  total  fractional  portion  of  the  sky  covered  by  clouds  (where  all  the  clouds  are  assumed 

*Aoo(a)  is  found  as  the  solution  to  the  infinite  scries  1  -f  a:  -f- =  (1  —  x)“\  where 
X  =  agQffff  is  that  portion  of  the  energy  trapped  between  the  foliage  and  surface  that 
survives  a  round  trip  of  reflection  at  the  surface  and  at  the  foliage  layer. 

^Since  these  emissions  only  approximate  a  perfect  blackbody  spectral  shape,  they  are  often 
called  graybody  radiation.  Gaseous  emitters  have  gaps  in  their  emission  spectra,  called 
windows,  where  no  emission  bands  exist.  Solid  surfaces  have  smoother  emission  spectra, 
but  emit  at  slightly  below  100  percent  efficiency. 


to  be  incorporated  into  a  single  parameter  at  a  single  hdght).  As  previously  discussed,  the 
temperature  of  the  upper  air  (cloud  temperature,  T,*)  is  used  to  characterize  the  longWTive 
emissions  source  region.  The  assignment  of  a  temperature  value  to  this  region  is  a  “best 
guess.”* 

The  a  used  in  equation  (17)  is  the  Stefan- Boltzmann  constant  (equal  to  5.67  x  10“®  W/(m^ 
K‘)). 

Note  that  the  cloud  temperature  is  represented  using  a  T  variable  rather  than  a  0  variable. 
The  T  variables  will  always  be  used  for  temperatures,  while  the  0  variables  will  always  be 
used  for  potential  temperatures.  The  two  are  related  through  the  equation 

0  =  T-Tz, 

where  z  is  the  height  above  the  ground  and  F  is  the  adiabatic  lapse  rate  (-9.8  K/km 
in  dry  atmospheres,  -6.5  K/km  in  moist  atmospheres).  For  a  neutrally  stable  (adiabatic) 
atmosphere,  0  would  be  constant  with  height  while  T  would  decrease  at  the  rate  of  F.  Since 
most  atmospheric  structure  models  are  dependent  on  the  stability,  d0/dz  is  much  more 
important  than  dT/dz,  and  thus  ^  is  a  more  natural  choice  for  the  temperature  variable 
near  the  surface  than  T.  But,  when  radiant  emissions  are  considered,  the  amount  of  radiant 
energy  depends  on  T  not  0.  Because  F  is  small,  assuming  measured  air  temperature  Ta 
is  approximately  equal  to  0a  represents  minimal  error  so  long  ^ls  the  station  height  z,  is 
within  a  few  meters  of  the  ground.  However,  for  the  longwave  emissions  from  clouds,  the 
parameter  T„a  must  be  used. 

The  atmospheric  emissivity  (ca)  used  in  equation  (17)  has  been  parameterized  by  Dcardorff 
(1978)  tising  the  theory  by  Staley  and  Jurica  (1972).  In  the  present  model  a  slightly 
modified  version  of  the  Deardorff  equation  is  used: 

Ca  =0.67e°  “®  (1.52 -0.00180a),  (18) 

where  c  is  the  atmospheric  vapor  pressure  at  the  station  height  in  millibars. 

The  deviation  of  equation  (18)  from  that  of  Deardorff  is  the  term  linear  in  0a,  which 
modifies  the  emissivity  based  on  temperature.  The  rationale  for  this  modification  is  that 
the  atmosphere  contains  two  window  regions  at  3-5  /im  and  8-12  ^m  in  the  infrared. 
Otherwise,  the  atmosphere  is  nearly  opaque.  If  one  therefore  assumes  that  the  atmosphere 
is  emitting  energy  as  a  blackbody  at  2J1  wavelengths  except  those  in  the  two  window  regions 
smd  that  no  emissions  in  the  window  regions  occur,  then  one  can  compute  the  fraction  of 
energy  emitted  (absorbed)  by  such  a  graybody  at  a  given  temperature.  Performing  this 
teisk,  one  finds  that  at  288  K  the  emissivity  of  such  a  graybody  is  0.74,  or  a  result  very 
similar  to  the  modeled  result  of  Staley  and  Jurica.  This  tends  to  indicate  the  physics 
of  the  approach  described  is  encompassing  the  essentials  of  the  emission  process  being 
modeled  by  the  original  equation.  Upon  performing  the  calculation  at  a  series  of  other 


*Al8o  notice  that  the  air  temperature  0*  used  in  the  longwave  flux  calculation  is  based 
on  the  input  air  temperature  rather  than  the  computed  temperature,  0,.  This  keeps  the 
downward  longwave  heat  flux  from  contributing  to  any  temperature  creep  effects. 
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temperatures,  one  obtains  a  curve  that  exhibits  a  nearly  linear  trend  with  temperatiure. 
This  trend  is  given  by  the  correction  term  in  equation  (17). 

Interestingly,  the  slope  of  the  correction  is  negative  since  with  increasing  temperature  a 
greater  fraction  of  the  radiant  energy  would  be  emitted  from  the  window  regions  (if  the 
windows  were  radiatively  active).  Since  the  windows  are  not  radiatively  active,  the  effective 
emissivity  decreases.  The  choice  of  coefficients  1.52  and  0.0018  is  such  that  the  correction 
term  equals  1  when  6a  =  288  K  (or  for  a  standard  atmospheric  temperature  similar  to 
that  tmder  which  the  results  of  Staley  and  Jurica  would  have  been  obtained).  At  higher 
temperatures,  eventually  the  effective  emissivity  rises  once  more,  but  the  temperature  at 
which  this  occvirs  is  far  above  those  normally  encountered  on  earth,  and  thus  the  linear 
trend  is  all  that  is  necessary. 

The  sky  emissions  have  therefore  been  divided  into  two  segments.  In  the  first  segment, 
the  near  surface  air  temperature  is  used  to  describe  emissions  from  air  near  the  surface  at 
wavelengths  outside  the  atmospheric  window  regions.  In  the  second  segment,  the  clouds 
are  treated  as  blackbody  surfaces  at  the  upper  air  temperature.  The  energy  arriving  from 
the  clouds  is  treated  as  proportional  to  the  fraction  of  the  sky  covered  by  clouds  and  as 
proportional  to  that  portion  of  the  energy  emitted  by  the  clouds  that  is  in  the  window 
region  of  the  spectrum  and  would  not  be  absorbed  by  the  intervening  air  on  its  way  to 
the  surface.  The  calculation  of  e^a  is  therefore  the  same  as  equation  (18),  except  that  Tua 
replaces  6a  since  the  spectrum  of  energy  depends  on  the  cloud  temperature  and  not  the 
near-surface  air  temperature. 


2.1.5  Other  Terrestrial  Radiative  Fluxes 

The  remaining  radiative  fluxes  originate  due  to  reflective  and  emissive  properties  of  the 
foliage  (/)  and  ground  {g)  surfaces.  To  represent  these  quantities,  Deardorff  used  Oj 
and  Of  for  the  mean  albedos  of  the  ground  and  foliage  surfaces,  respectively,  at  visible 
wavelengths.  Cj  and  tj  were  assigned  the  average  longwave  radiative  emissivities  for 
ground  and  foliage  (sec  Deardorff  (1978)  or  Oke  (1978)  for  tables  of  typical  values).  (Thus 
the  quantities  (1  —  Cj)  and  (1  —  €/)  represent  the  longwave  reflectivity  coefficients  for  the 
ground  and  foliage.) 

FVom  these  quantities  Deardorff  estimated  the  remaining  six  radiative  fluxes  in  equations 
(2)  and  (3).  Equations  (19)  and  (20)  are  similar  to  the  equations  Deardorff  used  for  the 
shortwave  and  longwave  fluxes. 


—  (1  •^oo(»)> 

(19a) 

5jt  =  OCg  Sg^  , 

(19i) 

Sh^  =  {1  ~  <Tf)Sg^  +  (T/af  Sh^  =  OeSki, 

(19c) 

Rlgi  =  (1  -  Aoo(t) 

(20o) 

+  <Tf  (1  -c/)e,a0j]Aoo(o> 

Rlgt  =  (1  -<rf)[Cg<T&*g+{l-eg)RLki]Aoo{t) 

(206) 
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+  a/  [c,  <r  ^  +  (1  -  e,)  tj  a  $)]  Aoo(t), 

Rhh^  =  (1  -  <rf)  ^  +  (1  -  ««)  [(1  -  <rj)  Riki  +<rf€f<7  0^f]}  Aoo(/)  (20c) 

+  {^/  ^  +  (1  “  */)  Rl^i  }» 

Aoo(i)  =  (1  -  (1  -  ef)(l  -  6,) a/)-».  (20d) 

Note  that  in  equation  (19a)  the  amovint  of  energy  reaching  the  surface  is  the  fraction  of 
direct  energy  passing  through  the  foliated  layer  multiplied  by  the  factor  >loo(«)-  The  factor 
■Aoo(*)  was  introduced  in  equation  (15b).  Likewise,  Aoe,(e)  represents  the  effects  of  multiple 
reflections  in  the  infrared.  The  main  assumption  made  Ixere  is  that  the  vegetation  layer  can 
be  effectively  replaced  by  a  partially  opaque  plane  surface.  Actually,  this  type  vegetation 
layer  will  never  exist,  but  mathematically  it  allows  for  a  simple  calculation  of  the  multiple 
reflections.  The  Deardorff  model  assumes  that  only  single  scattering  occurs.  Ekjuations 
(19)  and  (20)  are  thus  a  small  improvement  on  the  Deardorff  model.  Note,  however,  that 
normally  Aoo(a)  will  be  close  to  unity,  so  that  Deardorff ’s  original  equations  obtain  nearly 
the  same  results. 


2.1.6  Convective  Heat  and  Moisture  Fluxes 

The  fluxes  of  sensible  and  latent  heat  are  driven  by  the  bouyant  and  forced  mixing  of  a 
spectrum  of  different  sized  turbulent  air  elements  close  to  the  ground.  The  general  term 
of  “eddy”  is  used  to  describe  these  elements.  Many  simplified  models  have  been  used 
to  calculate  these  fluxes,  of  which  Deardorff’s  model  is  typical.  These  models  are  called 
simplified  because  a  fixed  value  is  used  for  the  heat  exchange  coefficient  during  all  unstable 
(daytime)  cases  and  a  second  fixed  quantity  for  all  stable  (nighttime)  cases.  However,  the 
Deardorff  model  also  considered  the  incorporation  of  a  foliated  layer,  and  this  aspect  of 
his  model  (suitably  modified)  has  been  transitioned  to  the  model  discussed  here. 

To  model  the  vegetation  layer,  Deardorff  used  a  modified  air  temperature,  vapor  pressure, 
and  windspeed  that  were  to  be  characteristic  of  the  values  of  these  parameters  in  the 
vicinity  of  the  leaves  of  the  vegetation.  While  some  other  aspects  of  the  Deardorff  approach 
will  be  retained,  this  particular  aspect  (of  using  near  leaf  parameters)  will  be  modified 
because  it  leads  to  overprediction  of  sensible  heat  for  highly  foliated  surfaces. 

A  second  (and  somewhat  better)  method  for  calculating  the  sensible  and  latent  heat 
fluxes  is  through  the  flux-profile  technique  employed  by  Hoffert  and  Storch  (1979)  for 
nonvegetated  surfaces.  Once  the  Hoffert  and  Storch  model  for  sensible  heat  flux  is  suitably 
modified  for  foliated  surface  effects,  constant  flux  coefficients  are  unnecessary  since  the 
method  develops  its  own  coefficient  values  dynamically,  based  on  the  current  atmospheric 
stability. 

The  basis  of  the  modified  approach  detailed  here  is  the  flux-profile  method.  The  flux- 
profile  method  uses  scaling  parameters  that  characterize  the  current  atmospheric  state 
within  the  surface  layer  (Dyer,  1974).  The  scaling  parameters  are  obtained  through  flux- 
profile  relationships.  These  parameters  can  then  be  related  to  vertical  profile  shapes,  based 
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on  the  similarity  theory  of  Obukhov  (1946).  Using  this  technique,  the  sensible  heat  flux 
(Hsh)  is  determined  by  the  equation 


Hsk  — —pCpT^Ut,  (21) 

where  p  is  the  dry  air  density,  Cp  is  the  specific  heat  of  dry  air  at  constant  pressure,  T»  is 
the  scaling  temperature  (Hoflert  and  Storch,  1979),  and  u,  is  the  friction  velocity. 

Using  this  equation  for  sensible  heat  flux  is  a  departure  from  the  DeardorfF  method,  which 
uses  an  equation  similar  to  that  outlined  by  Priestley  (1959,  see  pp.  4-6  and  39-43).* 

The  justification  behind  the  search  for  a  more  satisfying  calculation  for  sensible  heat  can 
be  traced  back  to  the  initial  purpose  of  the  model.  Originally,  the  author  was  looking  into 
the  atmospheric  refraction  problem.  In  this  problem  the  dominating  factor  influencing 
the  daytime  vertical  temperature  profile  is  the  sensible  heat  flux  Hsh-  Errors  in  the 
estimation  of  Hsh  would  therefore  immediately  be  translated  into  errors  in  estimation 
of  the  temperature  structtire.  Initial  model  experimentation  using  Deardorff’s  method 
showed  sensible  heat  fluxes  up  to  700  W/m*  for  heavily  fohated  conditions.  These  high 
values  were  considered  erroneous  and  found  to  be  the  result  of  using  an  inflated  value 
for  sensible  heat  flux  flow  coefficient,  Kut  when  <t/  was  large.  The  justification  used 
by  DeardorfF  to  explmn  the  increase  was  that  the  exposed  surface  area  for  sensible  heat 
flux  exchange  was  increasing  directly  as  the  fractional  foliage  cover  increased.  Thus  the 
reasoning  followed  that  higher  fluxes  would  result  from  a  higher  exposed  surface  area. 

The  problem  with  this  logic  is  that  the  sensible  heat,  according  to  flux-profile  theory, 
depends  only  on  the  shaiies  of  the  vertical  profiles  of  temperature  and  winclspeed,  not 
on  the  exposed  area  near  the  surface.  Flux-profile  theory  thus  predicts  that  there  is  a 
saturation  effect  on  the  carrying  capacity  of  the  air  that  is  passing  through  the  foliated 
layer  and  next  to  the  surface.  The  air  must  travel  close  enough  to  the  leaves,  stems, 
and  ground  surface  to  transfer  energy  and  then  must  be  replaced  by  new  air  from  above 
that  has  not  yet  been  influenced  by  the  surface  elements.  But  the  air  near  the  surface 
cannot  be  immediately  replaced  by  air  &om  above.  Thus  some  portions  of  the  plant  and 
soil  surfaces  will  always  be  in  contact  with  air  that  has  become  partially  adjusted  to  the 
warmer  (or  cooler)  temperature.  DeardorfF  tried  to  consider  this  aspect  by  parameterizing 
temperatures  very  close  to  the  surface  and  leaves,  but  his  attempt  was  only  partially 
successful.  In  fact  the  more  foliage  present,  the  more  of  a  barrier  will  exist  to  restrict  the 
free  flow  of  air  into  the  vegetated  canopy,  and  the  greater  effect  the  drag  of  the  foliage  will 
have  on  the  overall  windflow  patterns  near  the  surface. 

DrHi)it<r  this  limitation,  the  DeardorfF  model  documented  the  increased  efficiency  of  a 
foliat<'(l  lay<T  in  <;xchmiging  heat  and  moisture  as  compared  to  a  nonfoliated  surface. 
Allen  and  Lemon’s  (1972)  findings  were  cited  which  showed  that  for  each  square  meter 
of  overshadowing  foliage  there  will  be  about  7  m*  of  leaf  surface.  DeardorfF  complemented 


*DeardorfF  uses  an  approach  similar  to  Priestley’s,  except  that  the  coefficients  of  heat  and 
evaporation  (Priestley’s  Kh  and  Kw)  were  treated  as  functions  of  the  fractional  foliage 
cover  and  atmospheric  stability. 
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this  factor  of  7  by  an  additional  10  percent  to  account  for  “stalks,  stems,  twigs,  and  limbs 
which  exchange  heat  but  do  not  transpire.” 

In  the  Deardorff  approach  the  7.7  factor  scaled  directly  with  the  sensible  heat  flux.  Thus 
a  1  K  temperature  difference  between  the  ground/foliage  and  the  £iir  would  be  about  8.7 
(7.7  plus  the  1  for  the  groimd)  times  more  effective  in  transfering  heat  into  the  air  than 
would  a  bare  surface  with  a  similar  temperature  differential.  In  the  new  methodology,  a 
flux  profile  technique  is  modified  to  account  for  a  foliated  layer.  In  this  approach  Thom’s 
(1972)  findings  were  deemed  highly  applicable. 


&.1.7  Adapting  Thom’s  Resistance  Model  to  a  Calculation  of  Sensible  Heat 

In  his  paper  Thom  states  “The  aerodynamic  resistance  [rp]  encountered  at  a  rough 
surface  by  a  property  flux  will  in  general  exceed  the  resistance  encountered  there  by 
the  accompanying  flux  of  momentum  [ra].  This  excess  is  conveniently  expressed,  non- 
dimensionally  by  the  parameter  B~^  ...  for  the  property  ‘P’.” 

Bp^  =  u*(rp(z)  -  r«(2)).  (22) 


In  perhaps  more  plain  English,  Thom  is  speaking  about  a  technique  called  the  resistance 
methotl.  In  thf;  resistance  method,  the  property  fluxes  of  heat,  moment»jm,  and  humidity 
arc  treated  much  like  the  currents  (/)  in  electrical  circuits.  Similarly  the  temperature, 
humidity,  or  kinetic  energy  differences  between  the  air  at  station  height  and  at  the  surface 
can  be  likened  to  potential  differences  {V),  while  the  resistance  to  the  flow  (r)  is  a  parameter 
derived  from  surface  characteristics  and  the  wind  structure.  The  flow  rate  is  thus  found 
from  the  relation  I  =  Vfr. 

Thom  quantizes  the  bziseline  atmospheric  resistance  (resistance  to  momentum  flux)  as 


raiz) 


His  equation  for  rp  (the  environmental  resistance  to  the  flux  of  property  P)  is 


(23) 


rp(2) 


u* 


Xz  -  Xr 


X. 


(24) 


In  the  case  of  .sensibh*  heat,  Xt  is  the  temperature  at  height  z,  Xm  is  the  mean  surface 
temperature,  6^]*  and  x*  is  the  friction  concentration  of  temperature,  T*.  Equation  (24), 
when  applied  to  sensible  heat  fluxes,  is  thus  rewritten  as 


(25) 


*6e  8h2Jl  also  be  referred  to  as  the  effective  surface  temperature. 
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(The  iise  of  the  tenn  T*,  as  opposed  to  is  not  accidental.  It  has  historical  significance. 
T«  contains  a  k  within  its  definition  (see  Kunkel  and  Walters  (1983)  or  Hoffert  and  Storch 
(1979),  for  example),  while  does  not  (Paulson,  1970).) 

The  nondimensionalized  excess  resistance  to  heat  from  equation  (22))  could  now  be 
written  from  equations  (23)  and  (25).  However,  a  more  helpful  result  is  obtained  using  an 
additional  contribution  by  Thom.  For  Thom  developed  the  parameterization 

=  6.27  ui.  (26) 


Therefore,  combining  equations  (22),  (23),- and  (26),  r//  can  be  solved  directly. 


r/f(2) 


I 

«• 


(27) 


fVom  the  definition  of  T*,  the  ratio  within  the  absolute  value  sign  in  equation  (25)  will 
always  be  positive.  The  absolute  vaJue  may  therefore  be  removed,  and  T*  can  be  solved 
for  directly  from  equations  (27)  and  (25). 


(^z  ^e) 

(B^‘«,  +  u)’ 


(28) 


The  T*  thus  developed  is  applicable  to  a  fully  vegetated  layer. 

At  this  point  the  definition  for  T,  could  be  used  to  calculate  the  sensible  heat  in  equa¬ 
tion  (21).  However,  the  interactions  for  a  nonfoliated  surface  should  be  considered  and 
merged  with  result  (28)  for  foliated  surfaces  so  that  a  final  equation  will  apply  for  any  level 
of  foliation. 

In  this  regard,  the  results  of  Kunkel  and  Walters  (1983)  were  specifically  tmlored  for  a 
nonfoliated  surface:  Their  results  included  a  laminar  sublayer  pmameterization  to  handle 
a  completely  barren  surface. 

The  Kunkel  and  Walters  form  for  sensible  heat  flux  is  written  as 


Bsh  =  pCpUCH  (de  -  ^z)‘ 

By  comparison,  the  Thom  equivalent  form  for  sensible  heat  flux  is  written 


(21a) 


Bsfc  —  pCp 


{Bz-Bz) 


The  definition  for  6^  used  in  equation  (21a)  is 


*“  (l  +  l.lA^)  ’ 


(216) 


(29) 
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where 


N  =  7a  f.  (30) 

The  coefhcient  value  of  7  was  derived  by  Deardorff  from  findings  reported  by  Allen  and 
Lemon  (1972)  (for  a  com  crop)  and  by  Monteith  ci  al.  (1965)  (for  a  barley  crop)  that 
for  each  square  meter  of  overshadowing  foliage,  there  will  be  about  7  m^  of  leaf  surface 
area.  The  factor  1.1  represents  an  additional  10  percent  surface  area  that  accounts  for 
“stalks,  stems,  twigs,  and  limbs  which  exchange  heat  but  do  not  transpire.”  The  1.1 
factor  therefore  appezirs  in  the  sensible  heat  flux  equation,  but  not  in  the  latent  heat  flux 
equations. 

The  effective  temperature  of  the  surface  thus  represents  a  weighted  average  (based  on 
exposed  surface  area)  of  the  temperatures  of  the  soil  (through  0g)  and  the  foliage  (through 
0f).  Since  the  surface  area  of  the  leaves  can  be  proportionally  much  greater  than  that 
of  the  soil,  the  influence  of  the  groimd  will  be  considerably  reduced  when  a  heavily 
foliated  layer  is  present.  The  use  of  an  effective  surface  temperature  in  equation  (29) 
represents  the  means  used  to  correct  the  equation  for  heat  flux  proposed  by  Deardorff  for 
the  overestimation  problem  mentioned  previously.  The  weighting  function  allows  the  heat 
flux  to  be  proportionally  influenced  by  the  foliage,  but  the  overall  level  of  the  flux  will  not 
be  vastly  greater  simply  due  to  the  increased  surface  area  presented  by  the  foliage. 

Simple  comparison  of  equations  (21a)  and  (21b)  shows  the  only  difference  between  the 
Kunkcl/Waltcrs  and  Thom  equations  for  sensible  heat  are  the  c//  u  factor  used  by 
Kunkcl /Walters  and  the  factor  used  by  Thom.  Let  us  now  consider  in  some  detail  the 
functional  form  of  thcs<;  two  approaches  since  the  simple  appearance  of  c//  is  deceptive.  In 
Kunkel/Walters,  the  form  used  for  the  c//  equation  is 

"  0.74  (In  ^  -  V'm)(ln  f^-xkh^k  V fQ.l A)  ’ 

where  T'  was  derived  from  the  bluff  body  form  of  Garratt  and  Hicks  (1973),  V’m  ?ind  V’fc  are 
the  diabatic  influence  functions  described  in  section  2.4,  z  is  the  height  above  the  surface 
(which  in  this  case  will  be  set  to  the  station  height  z,),  Zq  :s  the  roughness  length,  and  k 
is  von  Karman’s  constant  (normally  set  to  0.4  for  rough  surfaces).* 

T'  is  given  by 

T'  =  [0.37(30u,zo/»')“‘‘®^r°®]. 

where  v  is  the  kinematic  viscosity  and  Pr  is  the  Prandtl  number  (equal  to  0.72  for  air). 
The  kinematic  viscosity  is  a  function  of  the  dynamic  viscosity  (ft)  and  the  density  of  air 
(p),  defined  as  (Smithsonian  Meteorological  Tables,  1951) 

u  =  n/p. 

*Kunkel  and  Walters  use  the  terms  ln(z/zo)  +  tfm  and  ln(z/zo)  +  ti>k^  but  neither  of  these 
usages  are  common,  and  in  their  equations  they  merely  use  deflnitions  for  xf;m  and  xph 
that  are  the  negative  of  those  commonly  used.  Therefore,  their  terminology'  has  been 
normalized  to  that  of  the  common  usages. 


22 


The  Smithsonian  Meteorological  Tables  also  define  n  as 


296.16  +  120  T®/* 

^  ”  r  +  120  296.16®/2  “  12.247  (T  + 120)  ’ 

where  T  is  the  air  temperature  in  Kelvin,  and  p  is 

p  =  0.34838  P/T, 

where  P  is  the  pressure  in  millibars,  and  the  resulting  density  is  in  kg/m®.  A  simple 
approximation  to  the  density  as  a  fimction  of  height  above  sea  level  (a.s.l.)  for  altitudes 
less  than  15  km  a.s.l.  has  been  empirically  derived  by  A.  Blanco  of  BE. 

,  =  1.225- ^[l.l76-j^[4.34  -  7.46j^]], 

where  H  is  the  height  a.s.l.  in  kilometers. 

The  equation  for  k  T'  therefore  reduces  to 

kT'  =  0.5258  (31) 

Since  the  air  referred  to  in  the  equations  that  require  v  is  only  a  few  millimeters  away 
from  the  ground,  6g  may  be  safely  substituted  for  T  in  the  density  and  dynamic  viscosity 
equations  above. 

Let  us  now  compare  the  terms 

C//U  ^  l/rniz), 

where  the  ^  symbol  is  used  to  denote  a  comparison  operation,  where  like  quantities  in 
two  dissimilar  developments  are  used  on  either  side  of  the  symbol.  If  the  Thom  form  for 
rfj  and  the  Kvmkel  and  Walters  form  for  ch  are  introduced,  we  have  the  relation 

O.74(ln^-0fc)-»-^r'  ^  u/u.  + V  ~ 
where  the  definition  (Paulson,  1970) 

u,  =  k  u/(ln(z/zo)  -  V»m) 

was  introduced  on  the  left-hand  side.  Now,  canceling  u,’s  on  both  sides  of  the  equation, 
dividing  by  k,  inverting  both  sides,  and  again  utilizing  the  definition  of  u*,  one  finds 

0.74(ln --tl^H)  +  kr  ^On-- 0m)  +  kBj,K 
Zo  Zo 

The  use  of  the  0.74  factor  was  once  believed  necessary  to  compare  the  parameterized 
temperature  structure  to  the  parameterized  wind  structure  (Businger,  1973),  but  this  was 


based  on  errors  committed  during  the  reduction  of  a  set  of  tower  data.  Upon  reexamination 
of  the  data  (Wieringa,  1980),  theorists  have  tended  to  set  the  0.74  constant  bsM;k  to  1.00, 
as  it  was  according  to  previous  theory.  Thus  the  log  terms  can  be  removed  and  the  result 
is 

Once  kT'  and  kBJi  are  written  as  functions  of  u«,  the  interrelationship  between  these 
approaches  becomes  even  clearer. 

ibT'  =  0.5258  6.27 kuV^  =  • 

Thus  an  analysis  of  the  two  approaches  finally  leads  to  comparable  terms  that  scale  as 
fractional  powers  of  the  friction  velocity  in  both  cases.  This  analysis  revesJs  that,  for 
what  would  appear  to  be  dissimilar  surfaces  (bare  ground  and  foliated  layer),  similarities 
still  exist  in  the  wind  interactions.  Perhaps  most  of  the  difference  in  exponents  can  be 
attributed  to  the  ability  of  plants  to  yield  to  the  wind  without  being  displaced,  while 
soil  may  be  physically  transported  once  its  static  friction  force  is  overcome,  but  normally 
remains  in  a  fixed  orientation  relative  to  the  wind  fiow. 

Regardless  of  these  differences,  the  results  of  the  two  approaches  are  similar  and  a 
combination  of  the  effects  of  both  responses  is  desired.  To  accomplish  this  combination  a 
mixing  coefficient  is  necessary.  A  linear  mixing  based  on  the  fractional  foliage  cover  was  at 
first  considered.  However,  since  the  roughness  length  zo  increases  rapidly  as  the  roughness 
elements  increase  in  size,  a  linear  mixing  model  would  apply  a  considerably  higher  surface 
roughness  length  to  the  k  T'  term  than  would  be  justified  by  the  conditions  under  which 
the  equation  itself  was  developed  (only  designed  for  smooth  surfaces  with  virtually  no 
foliage).  In  considering  the  form  for  the  mixing  rule  and  the  rule’s  relationship  to  the 
roughness  length  (which  depends  on  the  siuface’s  drag  effects  on  the  wind  profile),  a  new 
means  of  approximating  zq  was  also  considered. 

Ideally,  the  roughness  length  should  depend  on  the  type,  height,  and  fraction  of  the  surface 
covered  by  foliage  and  on  the  type  and  size  of  the  grains/i>ebbles/rocks  on  the  ground 
beneath  the  plants.  For  suitable  drag  effects,  the  following  rule  is  proposed  as  an  ad  hoc 
solution,  based  on  the  fractional  foliage  cover  and  estimates  of  the  roughness  lengths  due 
to  foliage  (zo/)  and  ground  below  (zoj)- 

Zo  w  Zo/  +  (1  -  <^y^)^og, 

where  zog  zq/.  For  example,  a  grass  covering  (roughness  elements  2”  to  5”  high  implies 
Zo/  «  1  cm)  over  a  sandy  soil  (roughness  elements  1/16”  to  1/2”  high  implies  zog  w 
0.06  cm)  with  o’/  =  0.5  would  yield  an  equivalent  roughness  length  of  0.81  cm.  In  essence 
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this  rule  states  that  the  roughness  parameter  rapidly  assumes  a  value  characteristic  of  the 
tallest  roughness  elements  as  the  fraction  of  those  dements  over  the  surface  increases.* 

Therefore  a  parameter  is  used  in  the  mixing  rtile,  rather  than  a/,  to  allow  for  a  much 
more  rapid  transition  from  a  smooth,  nonfoliated  surface  to  a  foliated  surface  dominated 
by  interactions  with  the  plants  rather  than  the  ground  surface. 

The  resultant  equation  for  T}i{z)  uses  a  linear  mixing  rule  with  a^'  as  the  mixing 
parameter  and  the  two  terms  —  ijim)  and  (kT*  —  0*)  as  the  terms  to  be  mixed. 

Once  mixed,  they  are  replaced  in  the  ori^nal  Thom  equation  for  resistance,  yielding 

•“(*)+  "y'  (*  +  (1  -  (h  r  -  V’.) 

nM  =  - - - r- - - - •  (32) 


Using  the  resistance  form  for  the  sensible  heat  flux  above  a  mixed  layer  (equation  (21b) 
with  equation  (32)  for  r//),  Hsh  can  be  effectively  divided  into  its  components  due  to  the 
surface  and  the  foliage,  respectively,  by  using  the  definition  of  6^ 


8,-e, 


(8.  -  e,) 

(l  +  l.lN) 


+ 


I.IN  \ 

1  +  i.iJvy  • 


Thus, 


Hsg  = 
Hsf  = 


Hsh  =  Hsg  +  Blsfy 

pCpi8,-8,)]  (  1  \ 

r„  JVl  +  l.lAr;’ 

>Cp(g/-g.)1  /  l.liV  \ 

r„  JV1  +  1.1N;’ 


In  the  program  these  relations  are  simplified  through  the  use  of  the  constants 


(33o) 

(336) 

(33c) 


Aa  = 


PCP  (  1  \ 

r„  Vl  +  l.lA^y’ 


(34a) 


A^^  =  l.lNAa,  (346) 

such  that 

Hsg  =  A»i8,-8,),  (35a) 


*The  accuracy  of  the  1/3  exponent  used  may  be  disputed,  but  the  point  is  clear:  In 
attempting  to  meld  the  two  approaches  it  is  necessary  to  avoid  applying  the  smooth  surface 
result  of  Kunkel  and  Walters  to  a  case  where  the  roughness  length  is  much  longer  than  the 
conditions  for  which  the  term  kT'  was  developed.  Otherwise  0.5258  (zo/i')®’^®  ^  6.27  k, 
and  the  .smooth  stirfacc  correction  term  may  exceed  the  foliage  correction  component, 
which  is  physically  unrealistic. 
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Hsf  =  AUBf-e,).  (356) 

In  a  similar  manner,  individual  latent  heat  fluxes  from  ground  and  foliage  are  found  in 
the  next  section.  Notice,  however,  that  the  sensible  heat  fluxes  from  soil  and  foliage  are 
assumed  to  be  completely  independent.  In  the  real  world  the  plants  could  possibly  heat 
up  first  and  thereby  warm  the  groxmd.  This  effect  is  not  modeled.  The  assumption  is  that 
any  interactions  between  groimd  and  foliage  will  be  through  radiant  fluxes  rather  than 
through  sensible  or  latent  heat  fluxes. 


2.1.8  Latent  Heat  Flux  for  a  Foliated  Layer 

In  the  Dcardorff  approach  to  moisture  transfer,  only  that  fraction  of  the  surface  that 
is  considered  to  be  moiat  may  be  involved  in  evapotranspiration.  The  parameters 
controlling  these  fractions  are  r"  and  o'  for  the  fractional  moist  foliage  and  ground  surfaces, 
respectively.  The  weighted  mean  fraction  of  the  surface  susceptible  to  evaporation  can  thus 
be  modeled  through  the  parameter 


o'  +  r"iV 
(l  +  N)- 

The  terms  o'  and  r"  are  based  on  Deardorff  (1978)  and  given  as 


o'  =  (1  -  6e)  +  min(l,  Wg/wk)6cy 


(36) 


r"  =  \~  6c  (r,/(r,  +  r^)]  [1  -  {tOdtw/wdmaxf^% 

^  _  f  1,  evaporation  is  occurring, 

~  (  0,  condensation  is  occurring, 

where  xug  is  the  fractional  moisture  of  the  ground  (by  volume);  lo/t  is  the  fraction  of  moisture 
the  ground  contains  when  it  behaves  as  if  it  is  saturated;  r,  is  the  stomat^Ll  resistance;  Tq 
is  the  atmospheric  resistance;  Wdew  is  the  mass  per  unit  ground  area  of  dew  on  the  foliage 
(normally  zero  during  the  day);  and  Wdmax  is  the  maximum  dew  accumulation  before  rtmoff 
to  the  ground  will  occur.  Further  equations  for  r,,  r^,  and  Wdmax  can  be  found  in  Deardorff 
(1978).  These  equations  depend  on  the  foliage  type,  wilt  factors,  mean  incident  sunlight, 
and  other  coupling  factors,  which  will  not  be  explained  here.  The  controlling  parameters 
are  generated  in  the  model,  and  no  justification  for  the  validity  of  Deardorff’s  approach 
will  be  presented.  Note  that  evaporation  is  considered  to  be  occurring  if  Q,at{Bg)  >  Qz  or 
>  Qzi  depending  on  whether  the  ground  or  foliage  surfaces  are  being  considered. 

Using  a  flux  equation  form  similar  to  that  in  Paulson  (1970),  the  total  latent  heat  flux  can 
be  described  by  the  equation 

L„Eh  =  -p  Lv  Q.  u*.  (37) 

Or,  using  the  same  method  used  to  develop  the  heat  equation, 


LvF^h  —  p  Li, 


{Q<-Qz) 


(38) 
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where  rq  =  rn,  since  all  fluxes  other  than  momentiun  flux  are  due  to  conservative  passive 
additive  transport  processes  and  the  same  dynamics  should  regulate  each.  Note  in  these 
flux  equations  that  the  A'  term  used  by  Deardorff  is  handled  implicitly  in  the  Qt  —  Qz 
difference  term. 

Having  deflned  A’  (equation  (36)),  and  since  Lv  is  the  latent  heat  of  vaporization  (equal  to 
2.5008  X 10*  J/kg)  and  Qt  is  the  speciflc  humidity  of  the  air  at  height  z  (Q,  =  Rh  Q»at{Tz) 
[see  appendix  A  for  computation  of  Q#ai(r)l),  *dl  that  remains  is  to  determine  Q*-  The 
used  in  equation  (38)  must  be  related  to  the  Q^at  values  of  the  leaf  and  ground  surfaces. 
Using  a  combination  of  equations  derived  by  Deardorff  and  the  same  weighted  averaging 
technique  used  in  the  sensible  heat  flux  derivation,  one  can  model  the  appropriate  surface 
effective  specific  humidity  as 


Qe  = 


a‘Qsaiie,)  +  r"NQta,{Bf) 

(l  +  JV) 


(39) 


This  weighting  method  emphasizes  that  the  equivalent  specific  humidity  of  the  surface  is 
related  to  both  the  relative  availability  of  water  on  the  ground  and  leaf  surfaces  and  the 
proportional  area  of  the  surface  occupied  by  the  given  surface  type.  Also,  the  1.1  term  is 
not  used  because  the  stems,  twigs,  and  so  forth  do  not  transpire.  The  weighting  achieved 
is  similar  to  weighting  the  Q,atiBx)  by  the  relative  humidity. 

With  these  definitions  for  the  quantities  govenung  the  overall  flux  of  moisture,  the  total 
latent  heat  flux  may  be  broken  into  its  respective  ground  and  foliage  related  components 


L^Eh  —  LtxiEg  “f*  LvEf^ 

(40a) 

where  the  constituent  fluxes  are 

^  ^  P —  QtIA')  / 

"  rq  \ 

.1  +  n)' 

(406) 

"  ^  rg  [ 

^  r"N  \ 

(40c) 

These  relationships  can  be  simplified  through  the  use  of  the  constants  Bt  and  B^ 

(41a) 

(416) 

so  that 

L„Eg  =  Bt{Q.,t{9t)-QzlA'l 
L„Ef  =  B^{Q,,,{ef)-QtlA'). 
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(42a) 

(426) 


2.2  Temperature  Calculations 

The  methods  of  evaluating  the  various  radiative  and  convective  heat  fluxes  were  discussed 
in  the  previous  section.  In  this  section  the  means  of  estimating  the  foliage  and  soil 
temperatures  are  described. 

The  Newton-Raphson  method  is  used  to  search  for  the  foliage  temperature  solution  to 


equation  (2),  which  can  be  rewritten  as 

Aal  5*,  +  Aa2  Rlhi  +  ^  —  -404  ^  =  Hsf  +  L„E/,  (43) 

where 

Aai  =  (1  -  at)  -  (1  -  as)(l  -  <rf)Aoo{i),  (44a) 

Aa2  =  1  +  <T/(1  -  C,)(l  -<T/)Aoo(|)  -«T/(1  “  «/)  “  (1  “  ^/)  ^oo(/),  (446) 

Aaz  =<Teg<Tf  tf  Aoo(0,  (44c) 

i4a4  =  <T  e f  (Tf  [I  +  Aoo(/)  -<Tf(l~  €g)  Aoo(i)l.  (44d) 


To  elucidate  further  on  the  components  of  the  Aai  constants:  Aai  has  two  terms.  The 
first  term  represents  the  incident  direct  solar  (the  1)  and  the  reflected  direct  solar  (— a;«). 
The  second  term  represents  the  energy  that  passes  through  the  layer,  —(1  —  (Tf),  and  the 
energy  that  returns  after  reflection  from  the  soil,  Og  (1  —  <t/). 

Aa2  has  four  terms.  The  first  term  represents  the  incident  longwave  radiation  from  the 
sky.  The  second  term  represents  the  energy  travelling  upward  toward  the  foliage  layer 
following  an  initial  pa.ssagc  through  the  layer  and  a  reflection  at  the  surface.  The  third 
term  represents  th<'  (uiergy  that  travels  upward  from  the  top  of  the  foliage.  The  fotjrth 
term  represents  the  energy  that  travels  downward  from  the  bottom  of  the  foliage  layer 
toward  the  surface. 

Aaz  hsis  only  the  single  term  related  to  absorption  by  the  foliage  layer  of  energy  emitted 
by  the  soil. 

Aa4  has  three  terms.  The  first  term  relates  to  the  upward  thermal  emissions  from  the 
foliage.  The  second  and  third  terms  refer  to  the  energies  that  flow  between  the  foliage 
layer  and  the  soil  that  originated  as  graybody  radiation  at  the  foliage  layer. 

If  the  constant  terms  are  gathered  together,  a  single  function  of  0f  may  be  formed.  Let 

Ail  =  Aai  Sh^  +  Aa2  Rlki  +  AaZ  0^  +  Aaa  0t  +  Bbb  QxiA' .  (44c) 

Equation  (43)  can  then  be  written  as 

F{0f)  =  -Ail  +  >l«4  ^  +  Aaa  0f  +  Bii  QaatiBf)  =  0.  (45) 

To  solve  for  F{0f)  =  0  using  Newton-Raphson,  the  previous  value  of  foliage  temperature 
is  used  as  the  initial  value  in  the  iterative  equation 


^/(n+I)  —  ^/(n) 


L 

F' 


(46) 
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until  the  equation  converges,  where  the  term  F'  is  (pven  by 


+  (47) 

The  definition  of  Q,at  and  its  derivative  are  discussed  in  appendix  A. 

In  this  derivation  the  assumption  is  that  6^  is  known  while  making  the  computations  to 
estimate  6f.  But,  we  cannot  know  dg  unless  we  know  6j^  which  leads  to  a  circular  argument. 
This  dilemma  is  ‘‘solved”  in  this  model  by  using  the  9f  from  the  previous  time  step.  Such 
an  assumption  avoids  certain  problems  due  to  mathematical  instabilities  in  the  model  that 
might  otherwise  arise.  It  should  be  a  good  assumption,  as  long  as  the  time  increment  used 
in  updating  the  soil  temperature  is  reasonably  short,  since  the  soil  temperature  (which 
includes  an  inertia  term)  will  vary  more  slowly  than  the  foliage  temperatmre.  The  foliage 
temperature  calculation  then  has  a  imique  solution.* 

Once  the  foliage  temperature  has  been  estimated,  the  new  groimd  temperature  is  found 
using  the  Crank-Nicholson  method.  This  method  assumes  the  time  step  (At)  is  small 
enough  that  fluxes  at  time  cein  be  determined  based  on  the  temperature  value  at 

time  temperature  difference  —  ^j(„))  across  the  time  step.  For 

example,  the  flux  of  thermal  rsuliative  emissions  from  the  ground  is  proportional  to  the 
temperature  to  the  fourth  power.  This  is  approximated  at  the  new  time  by 


^(»  +  0  ^  ^(»)  ^^(n) 


(48) 


Similarly,  the  new  saturated  specific  hvunidity  is  approximated  based  on  its  present  value 
and  a  change  with  temperature  of 


'  •»(") 

The  fluxes  during  the  period  covered  by  the  time  step  are  treated  using  average  values  for 
the  fluxes.  The  average<l  temperature  dependencies  for  sensible,  longwave,  and  evaporative 
heat  fluxes  arc 


9g  « 


^  ^  ^^(n)  ^9(n+l)  ^(»)» 


(50a) 

(506) 

(50c) 


*However,  when  reevaluating  F  and  F'  at  each  step  in  the  Newton-Raphson  approach,  it 
is  necessary  to  consider  changes  in  r",  since  when  on  the  borderline  between  evaporation 
and  condensation,  r"  can  vary  widely  depending  on  the  value  of  9/. 
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All  the  heat  fluxes  in  equation  (3)  can  now  be  written  as  functions  of  and 
These  quantities  may  then  be  substituted  for  G  in  equation  (4),  and  the  remaining  terms  in 
equation  (4)  can  also  be  expressed  as  functions  of  and  The  resulting  equation 

can  then  be  separated  into  that  portion  that  depends  on  and  on  the  remaining 

terms  that  are  considered  constant.  Since  the  resulting  function  is  linear  in  ^j(„+i) 

is  found  to  equal  the  ratio  of  two  sets  of  terms.  Following  this  approach,  one  c?in  rewrite 
equation  (4)  as 


(^9(n  +  l)  ^S(n))  _  Cl  G  ^  C2  ^  C2  +  ^g(n)  ) 


e, 


«("+») 


C2 


At  2Tday\ 


CiG  .  C2(2^2  .  ^J(„) 

I  ^  + 


h 


2Td, 


ay 


At  ’ 


(51a) 

(516) 


where  the  G  term  must  still  be  divided  between  constants  and  terms  linear  in  Bg^^^^y  This 
tjusk  is  accomplished  in  the  equation 


G  =  S„et  +RLg2  ~G$,  +  ^i 


where 


Snet  —  (1  ~  C(j)(l  —  <7/)  Shi  "d-ooC*)? 

Rlgi  =  tg  >loo(«)  {(1  -  <Tf)RLhi  +<Tf€f<r0)-[<Tfil-€f)-  }, 

C.„.,  =  (*,<., /2  -  «.)  +  B.  (^),_  -  ' 

»(n) 

Thus  the  solution  for  the  ground  temperature  value  at  the  next  time  step  is 


0 


CliS„el  +  RLg2  ^2  (2^2  ”  ^j(„) ) 

/l  2Tday  At 


9(n  +  l) 


1  C2  _  Cl  Fi 
At  2Tday  h  J 


(52) 

(53a) 

(536) 

(53c) 

(53d) 


(54) 


2.3  Modeling  Additional  Air,  Soil,  and  Plant  Parameters 

In  this  section  the  equations  related  to  92,  conservation  of  moisture,  and  are  discussed. 
Also,  parameters  that  evolve  with  time  due  to  their  dependence  on  moisture  are  considered. 
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i.S.l  MoisiuTt  and  Temperature  Equations 


In  addition  to  the  surface  and  foliage  temperatures,  the  deep  ground  temperature  62  and 
the  flux  dependent  near  surface  air  temperature  6^  must  be  computed.  The  first  of  these 
equations  is  quite  simple,  being  derived  from  equation  (6),  using  the  stepping  methodology 
devised  for  the  ground  temperature  that  was  discussed  in  the  previous  section: 

~  ^2(«)  _  £ 

At  I2 


whence 


B2 


"+» 


+ 


AtG 

I2  ■ 


The  net  flux  into  the  ground,  G,  is  known  as  a  by-product  of  computing  the  new 
temperature  (equation  (52)). 

A  similar  equation  to  that  developed  for  6g  is  used  to  determine  the  updated  values  for 
the  soil  moisture  parameters  Wg  and  W2.  The  change  in  Wg  is  modeled  as 


^9(n^X) 


—  w 


At 


^  =  -Cl 


(Eg  +  0.1  Ef-Precip)  ^^(ttfg-U>2) 


Pwd'j 


0<Wg<  tVjnaz,  (55) 


where  Predp  is  the  amount  of  precipitation  measured  in  kg/(m^  s),  dj  is  the  distance 
10  cm,  Wmax  is  the  maximum  fraction  of  moisture  the  near  surface  ground  can  contain 
before  runoff  occurs,  pw  is  the  density  of  water  (1000  kg/m®),  and  the  constants  Ci  and 
C2  are  weighting  terms  based  on  the  work  of  Jackson  (1973)  (see  Deardorff,  1978).  As 
seen  from  the  above  equation,  the  modelctl  surface  moisture  depends  on  transpiration, 
evaporation,  and  precipitation  in  the  first  term  on  the  right  and  on  a  restoring  flow  of 
moisture  up  through  the  ground  from  the  deeper  layer. 

Similar  to  the  computation  for  62,  the  computation  of  W2  depends  on  the  overall  flux  of 
moisture  through  the  surface. 


+  ^^(n)  _  (^9  -h  Ef  Prtcip) 

Att  Pw  d^2 

where  dj  is  the  distance  50  cm.  Notice  that  the  Wg  equation  only  incorporates  1/lOth 
of  the  Ef  transpiration  flux.  This  is  because  the  Deardorff  model  asstunes  the  roots  of 
the  plants  extend  deeper  than  the  surface  layer.  Note  also  that  the  constants  Ci  and  C2 
depend  on  empirical  rules  simulating  the  capillarity  and  other  transport  properties  of  the 
soil.  Deardorff  provides  no  information  other  than  single  constants  for  these  terms,  which 
have  no  direct  relation  to  physically  measureable  quantities. 

The  estimation  of  the  near  surface  air  temperature,  is  somewhat  more  problematic 
since  the  temperature  will  depend  not  only  on  the  measured  air  temperature  and  the 
estimated  sensible  heat  flux  but  also  on  the  characteristic  depth  of  mixing  expected  under 
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the  present  condition  of  the  atmospheric  stability. Based  on  these  assumptions,  the  mixing 
rule  developed  is 


e 


*(»+») 


Hsk^t 

pCpZinv' 


(57) 


where, 


=50u 


(58) 


is  an  ad  hoc  estimate  of  the  mixing  depth  of  the  atmosphere.  The  constant  720  reflects 
the  idea  that  the  maximum  lag  time  allowed  in  the  model  between  measured  temperature 
Ba  and  modeled  temperature  Bx  is  only  half  a  day  of  model  time  and  allows  for  seasonal 
and  synoptic  influences  on  the  energy  budget.  The  term  varying  with  sensible  heat  means 
the  near  surface  air  temperature  will  respond  to  the  diurnal  cycles  of  warming  and  cooling 
that  occur. 


2.3.2  Moisture  Dependent  Coefficients 


Perhaps  after  computing  the  fluxes  and  temperatures,  one  might  imagine  they  could 
completely  characterize  the  turbvilent  structures  and  vertical  mean  profiles.  This  is 
possible,  if  the  total  simulated  time  is  short.  However,  first,  the  model  requires  at  least 
a  few  hours  of  simulated  running  time  to  settle  down  from  its  initialization  (due  to  the 
thermal  inertia  term).  Second,  if  the  model  is  nm  so  as  to  simulate  a  few  days  time,  the 
soil  moisture  parameters  will  alter  the  soil  thermal  characteristics.  This  bookkeeping  can 
be  handled  between  the  temperattire  computations  at  each  time  step. 

Deardorff  gave  parameterized  equations  to  allow  the  heat  conduction,  heat  capacity,  and 
surface  reflectivity  to  vary  with  the  soil  moisture.  In  suggesting  a  modification  on  the 
Deardorff  approach,  one  might  allow  a  user  to  designate  a  set  of  standard  soil  types  or 
input  unique  initial  conditions.  As  the  model  runs,  these  variables  are  modified  by  the 
change  in  water  content  from  the  initial  conditions  as  follows: 


P»  c,  =  {p,  c,)o  +  4.184  X  10*  (u;  -  u’o), 


(59a) 


Pa  Cg  Ka  =  X  =  Xq  +  K,  X  4.184  X  10* 


o«  =  tto  —  0.17 


{Wg  -  WyJ 


Wk 


(596) 

(59c) 


where  p,  c,  is  the  soil  heat  capacity  (soil  density  times  the  specific  heat  capacity),  p,  c,  k, 
is  the  thermal  conductivity  (with  k,  the  thermal  diifusivity),  and  Qg  has  already  been 


*Note  that  latent  heat  is  not  expected  to  contribute  to  the  temperature  change  in 
the  atmosphere,  and  in  local  thermodynamic  equilibrium  there  should  be  no  significant 
radiative  contributions  cither. 


32 


given  as  the  soil  mean  reflectivity  in  the  shortwave  band  (albedo).  The  first  two  equations 
have  been  applied  to  the  surface  (wg)  and  deep  soil  (W2)  thermal  parameters.  These 
equations  parallel  DeardorfF’s  equations  (37),  (38a),  and  (40).  Different  standard  soil 
types,  with  associated  nominal  properties,  were  included  in  the  TGRAD  model  as  options. 
The  soil  properties  were  derived  from  average  soil  conditions  as  documented  by  Oke  (1978), 
Deardorff  (1978),  and  Link  (1979).  The  final  equation  is  applied  at  the  surface,  and  TGRAD 
utilizes  reflectivity  coefiBcients  derived  from  the  same  sources. 


£.S.S  Possible  Improvements 

In  addition  to  variation  under  changing  moisture  in  the  soil,  it  is  possible  to  configure  a 
program  to  treat  snow  cover  and  freezing  conditions.  Such  a  model  would  be  beneficial 
in  approximating  a  European  climate.  However,  the  scope  of  such  a  model  is  very  broad, 
covering  the  various  aspects  of  adding  a  new  layer  in  the  temperature  flux  computations 
for  the  snow,  of  adding  the  additional  considerations  of  conversion  of  water  to  ice  and 
ice  to  water,  and  of  adding  various  other  complexities.  As  such,  a  model  overing  these 
topics  would  be  a  subject  unto  itself.  The  model  described  here  is  not  designed  to  make  a 
transition  through  freezing  or  thawing  conditions. 

Note  also  that  the  temperature  of  precipitation  may  influence  the  soil  temperature  calcu¬ 
lation,  but  was  excluded.  From  a  practical  standpoint,  the  structure  of  the  atmosphere 
will  normally  be  neutral  during  any  precipitation  episode. 

Another  consideration  is  in  regard  to  the  seasonal  variation  of  the  stomatal  resistance  of 
I)lants  to  evai)otrauspiration.  For  the  stomatal  resistance  equation,  Deardorff  indicates 
that  “in  temperate  latitudes,  S  (the  stomatal  resistance’s  seasonal  dependence  factor)  is 
set  to  zero  during  the  growing  season  and  to  a  value  much  larger  than  unity  during  the 
rest  of  the  ye^lr.”  In  the  program,  S  is  set  to  50  except  during  a  European  growing 
season.  For  desert  conditions  S  is  assumed  to  remain  high  year-round,  since,  for  most 
plants,  survival  requires  high  transpiration  resistance.  Relatively  little  is  known  about  the 
stomatal  resistance  of  the  various  weeds  and  scrub  brush  normally  encountered  in  actual 
tactical  terrains  and  during  field  measurements,  since  most  studies  have  focused  on  crop 
characteristics.  A  comprehensive  study  of  the  properties  of  general  vegetation  would  be  in 
order,  to  truly  obtain  a  more  rigorous  model. 

The  present  model  is  designed  to  include  sequential  inputs  of  data  and  compute  the 
resulting  atmospheric  states  as  a  time  series.  However,  the  model  could  be  reconfigured 
relatively  simply  so  that  it  could  run  a  single  day  case,  reinitializing  the  soil  parameters  at 
an  appropriate  time  each  day.  This  step  is  possible  using  simple  bookkeeping  techniqiies, 
whose  details  need  not  be  considered  further  here.  The  advantage  of  such  a  modification 
would  be  to  allow  for  the  settling  down  of  any  transient  behavior  to  obtain  a  repetitive 
(cyclical)  response  of  the  model  to  a  single  day’s  input  data  set. 
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2.4  Bulk  Aerodynamic  Parameters 


In  the  derivation  of  the  equation  for  sensible  heat,  the  terms  u,,  V’m,  and  were 
introduced,  but  not  defined  mathematically.  The  definitions  are  provided  in  this  section, 
finally  completing  the  description  of  the  means  of  estimating  the  surface  heat  fluxes.  In 
addition,  the  purpose  of  the  model  was  not  merely  to  compute  the  ground  and  foliage 
temperatures.  Rather,  the  computations  were  a  step  toward  characterizing  the  vertical 
structures  of  temperature,  refractive  index,  refractive  index  structure  parameter,  inner 
and  outer  scales,  and  windspeed.  Therefore,  in  this  section,  the  flux  profile  forms  for  the 
vertical  structures  are  introduced  and  the  means  of  estimating  u*,  V’m,  and  are  given, 
along  with  the  vertical  structure  equations  for  the  other  parameters. 

One  of  the  principle  scaling  parameters  used  in  flux  profile  theory  is  the  Obukhov  length. 
To  estimate  the  friction  velocity  u,  and  aerodynamic  resistance  r//,  the  Obukhov  length 
must  also  be  computed.  This  length  is  defined  as 

_  pCpOxul 
kgHsk  ’ 


where  g  is  the  gravitational  acceleration.  The  friction  velocity  is  found  as  a  function  of 
the  Obukhov  length. 


ku 

(^n(j)  -V>m) 


(Gl) 


where  zo  is  the  roughness  height,  z  is  the  height  above  the  surface  where  u  is  to  be 
measured,  and  V’m  is  the  diabatic  influence  function  for  momentum.  Of  course,  once 
«*  has  been  estimated,  the  vertical  structure  of  windspeed  can  be  determined  simply  by 
inverting  equation  (61)  to  compute 


(62) 


The  correct  definition  of  z  must  be  handled  carefully  since  it  is  defined  with  respect  to  the 
vertical  structure  of  the  wind  profile  and  not  simply  with  respect  to  the  height  above  the 
surface.  Foliated  layers  introduce  a  correction  to  the  value  of  z  compared  to  the  actual 
height  above  a  physical  surface.  In  particular,  let  r'  be  the  actual  height  of  the  windspeed 
measuring  station  above  the  ground,  z,  is  related  to  z'  through  the  relation 

z,  =  -  D,  (63) 

where  D  is  called  the  displacement  height.  This  quantity  is  found  using 

D  =  0.70(zo/0.13).  (64) 
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This  equation  is  derived  as  follows:  D  is  modeled  as  being  approximately  7/lOths  of  the 
height  of  the  vegetation.  The  height  of  the  vegetation  is  approximately  «o/0.13.*  Since  we 
previously  redefined  the  roughness  length  for  partially  vegetated  layers,  the  D  parameter 
should  vary  smoothly  toward  a  smaller  value  as  the  fractional  foliage  cover  is  reduced. 
Figure  2  shows  a  representative  vertical  wind  structure  pattern  when  the  displacement 
height  is  incorporated  (after  Oke  (1978)). 

In  some  models  the  height  of  the  extrapolated  surface  value  for  temperature  is  also 
computed.  In  the  approach  used  in  this  model,  this  height  difference  is  accounted  for 
differently  through  the  difference  in  resistances  between  momentum  and  sensible  heat 
fluxes. 


Figure  2. 


Foliage  height  H  and  extrapolated  zero  of  the  windspeed  profile  for  a  fully 
foliated  layer. 


Returning  now  to  the  main  line  of  argument,  in  equation  (61)  the  term  V’m  is  the  diabatic 
influence  function.  (The  m  stands  for  momentum  and  is  usually  applied  as  a  modifier  to 
the  wind  profile.)  The  symbol  V’fc  similarly  represents  the  diabatic  influence  function  for 
sensible  heat  (h).  As  helps  scale  the  vertical  wind  profile,  V’A  helps  scale  the  vertical 
temperature  profile. 


V’m  and  V’A  are  both  functions  of  height  {z)  and  Obukhov  length  and  have  different 
mathematical  forms  depending  on  the  stability  conditions.  When  the  Obukhov  length 
is  less  than  zefo  (an  unstable  atmosphere)  the  V*  functions  have  the  form 


=  2£n 


—  2  tan  * 


(65a) 


*Frank  Hansen,  c.  1990,  BE,  WSMR,  NM,  informal  commimication,  indicates  the  rough¬ 
ness  length  is  usually  about  13  percent  of  the  height  of  the  roughness  elements.  Thus  by 
dividing  zo  by  0.13,  the  height  of  the  roughness  elements  is  obtained. 
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(655) 

=  2k)' 

(65c) 

(65d) 

When  the  Obukhov  length  is  greater  than  zero  (a  stable  atmosphere)  the  functions  take 
the  form 

il>m  =  =  I  -  <f>,  {06a) 

=  {OOh) 


The  terms  stable  and  unstable  refer  to  the  behavior  of  a  parcel  of  air  at  the  same 
temperature  with  height  as  its  surroundings  when  that  parcel  is  displaced  slightly  in  a 
vertical  direction.  In  a  stable  atmosphere  the  parcel  will  attempt  to  return  to  the  level 
from  which  it  was  displaced  (perhaps  initiating  an  oscillatory  motion).  In  an  unstable 
atmosphere  the  parcel  will  become  warmer  or  cooler  than  its  environment,  depending  on 
whether  the  displacement  was  upward  or  downward,  respectively.  In  this  case  the  force  on 
the  parcel  will  be  such  that  the  parcel  will  move  vertically  away  from  its  level  of  origin. 

Since  the  Obukhov  length  is  negative  for  unstable  atmospheres  and  positive  for  stable 
atmospheres,  it  is  a  useful  measure  of  the  current  state  of  the  vertical  structure  of  the 
atmosphere.  It  approaches  either  positive  or  negative  infinity  under  neutral  (adiabatic) 
stability  conditions.  Normally,  the  Obukhov  length  appears  in  equations  in  the  denomina¬ 
tor  of  a  dimen.sionless  length  variable  such  as  ^  =  z/L.  Thus  there  is  a  smooth  transition  of 
the  (  variable  across  the  neutral  stability  condition  between  stable  (positive)  and  unstable 
(negative)  atmospheres  as  L  varies. 

Although  there  are  several  other  forms  for  equations  (65)  and  (66)  available  in  the 
literature  (Yaglom,  1977),  the  forms  chosen  were  those  suggested  by  Hansen  in  pri\'ate 
communication*  and  correspond  to  the  equations  recommended  by  Dyer  (1974)  and  Hanna 
et  al  (1982). 

Notice  that  we  now  have  a  problem:  L  depends  on  u,  and  Hsk]  depends  on  L;  and 
Hsk  depends  on  «*,  r//,  df,  and  $g.  depends  (ultimately)  on  L  and  «*;  and  0/  and  0g 
depend  on  the  Hsk-  One  can  only  compute  ail  these  variables  simultaneously,  by  using  an 
iterative  procedure.  This  procedure  fixes  the  temperatme  difference  de  —  Bx  and  the  wind 
speed  u  to  constant  values  while  L,  r//,  ffsfc,  and  u*  are  successively  calculated  until  the 
Obukhov  length  converges.  The  foliage  and  surface  temperatures  can  then  be  computed. 
An  approximation  is  therefore  made  that  the  constants  Aaa,  Bt,  and  Bu  will  depend 
on  the  previous  time  interval’s  0g  and  0f  estimates.  As  long  as  the  time  interval  is  short. 


*The  ^  coefficient  is  equal  to  the  inverse  of  the  critical  Richardson  number. 
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this  approximation  should  not  lead  to  significant  error.  Certainly  it  is  less  in  error  than 
the  approximation  of  constant  flux  controlling  coefficients. 

To  simultaneously  calculate  the  above  parameters,  first  is  estimated  using  an  assumption 
of  an  adiabatic  atmosphere.  Then,  L  can  be  estimated.  Subsequently,  th  can  be 
estimated.  As  discussed  in  the  previous  paragraph,  6g  and  6/  will  be  used  in  determining 
the  temperature  difference  between  the  station  height  and  the  surface.  Hence,  the  estimate 
of  the  sensible  heat  flux  can  be  determined  immediately  after  estimating  r/f.  Once  these 
four  steps  are  complete,  the  cycle  can  be  repeated  to  determine  an  estimate  of  u*  with 
the  previous  (nonadiabatic)  value  for  L  and  more  accurate  values  of  the  other  parameters. 
The  cycle  is  halted  once  the  parameters  sufficiently  approach  limiting  values. 

In  following  this  procedure,  unstable  atmospheric  conditions  do  not  present  a  problem, 
and  computations  always  converge.  However,  for  stable  atmospheres  the  equations  may 
become  inconsistent.  For  the  moment,  ignoring  the  k  and  kT‘  influences  on  r//,  and 
since 

-rpm  =  -V'fc  =  <j>-l  =  ^  zIL, 

TH  can  be  approximated  by  u/uj-  The  equation  for  L  thus  reduces  to 

r  _  _ _ 

"  g{6,-6,){in{0  +  fizlLy 
where  —  Let  us  now  rewrite  this  equation, 

1  =  g-^MC)(^z  -ge)  ^gzje^-e^) 

6,  u2  ^ 

Notice  that  imder  stable  atmospheric  conditions,  L  must  always  be  greater  than  zero. 
Similarly,  $z  will  be  wanner  than  the  surface  temperature,  and  y,  zuid  z  are  all  positive 
values.  Therefore,  there  are  no  negative  terms  on  the  right  side  of  the  equation,  and  if 


I3gz{ez  -Br)  . 
Bzu^ 


the  equation  has  no  solution.  Thus  for  a  real  solution  the  condition 


g-l  9  zifit  —  Be) 

must  be  met.  When  this  condition  fails,  is  set  to  a  minimum  value  of  0.07  to  force 
a  reasonable  value  for  L  to  be  obtained.  With  u«  fixed,  r//  and  L  are  determined 
directly.  This  technique  is,  however,  less  than  optimal,  and  one  would  hope  for  a  better 
imderstanding  of  the  behavior  of  the  nocturnal  surface  layer  atmosphere.  Appendix  B 
includes  a  small  description  of  the  problems  encountered  at  night  and  some  anidysis  of 
gravity  wave  influence  on  the  overall  time-based  structure  of  the  nocturnal  temperature 
gradient. 
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The  CM  subroutine  also  computes  the  quantities  related  to  the  characterization  of  the 
surface  layer  atmosphere.  These  include  C*,  Cy,  Cy,  and  dTfdz.  Kunkel  and  Walters 
(1983)  provide  equations  for  the  C\  and  parameters. 

Cl  =  {C\.{z)A^  +  0.03/B)2,  (68) 

where  A  is  a  proportionality  constant  whose  value  depends  on  wavelength.  In  the 
visible/infrared  region  of  the  spectrum,  the  value  A  =  79  x  10“®K  mbar~^  is  applicable. 
The  quantity  B  is  called  the  Bowen  ratio  and  is  equal  to  HskILEh-  C^{z)  was  obtained 
from  Wyngaard  (1973)  as 


=  Tl  2“*/’  4.9  (1  -  70"*/’  (69) 

for  unstable  conditions,  and 

Cl  =  Tl  2"*/*  4.9(1  +  2.4^/’)  (70) 


for  stable  conditions. 

For  Cy,  Ochs  and  Hill  (1985)  recommend  the  equation 

Cl  ^  2€2/^  (71) 

where  the  Cy  computed  is  the  longitudinal  component  (the  most  significant  component 
for  acoustic  propagation  problems),  and  e  is  the  energy  dissipation  rate.  Businger  (1973) 
provides  an  equation  for  the  inner  scale  as  a  function  of  u,  and  z  under  neutral  stability 
conditions. 


Analysis  of  data  measured  by  Ochs  ttmds  to  indicate  the  dependence  of  e  on  stability 
conditions  is  weeik  (Tofsted  and  Auvermann,  1991).*  This  finding  would  appear  to  be  in 
accord  with  the  relationship  of  Cy  with  e  (Cy  is  most  dependent  on  the  turbulent  loss  of 
fluctuation  energy  due  to  viscous  dissipation.  But  viscous  dissipation  is  a  process  related 
to  kinetic  collisions  of  molecules  and  has  little  dependence  on  stability  at  this  small  size 
scale). 


*The  data  mentioned  was  unpublished  and  was  obtained  from  Gerald  Ochs  of  the  Wave 
Propagation  Laboratory  (WPL),  National  Oceanographic  and  Atmospheric  Administra¬ 
tion  (NOAA),  Boulder,  CO,  c.  1990.  The  data  was  a  windspeed  versus  inner  scale  scatter 
plot,  which  has  been  incorporated  into  Tofsted  and  Auvermann  (1991).  The  statement 
that  the  data  infers  minimal  dependence  of  the  dissipation  rate  on  stability  follows  from 
the  derivation  of  inner  scale,  which  depends  on  dissipation.  Had  there  been  a  marked 
dependence  of  dissipation  on  stability,  the  plot  would  have  had  significantly  greater  scatter 
than  was  present,  particularly  at  low  windspeeds. 
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For  acoustic  and  visible  refraction,  the  temperature  gradient  is  given  by 


where  F  is  -0.0098  K/m  for  relative  humidities  less  than  about  98  percent  and  -0.0065  K/m 
for  higher  relative  humidities. 

The  equation  for  the  turbulent  outer  scale  was  derived  from  a  combination  of  the  param¬ 
eterization  for  dufdz,  obtained  from  Paulson  (1970),  and  an  equation  for  outer  scale  that 
depends  on  du/dz  obtained  from  LeweUen  (1977). 

Io  =  1.68  ibz  (1-16  C)‘^^  (74) 


where  k  is  again  von  Karmw’s  constant  and  ^  =  z/I». 

Finally,  an  expression  for  inner  scale  was  derived  in  Tofsted  and  Auvermann  (1991).  This 
development  took  an  analytical  expression  for  inner  scale  (Wyngaard,  1973)  and  modified 
it  by  an  empirical  expression  based  on  data  measured  at  Table  Mountain,  Colorado,  by 
Ochs.  The  expression  was 


£o  «  0.000463 


z^/^en(z/zQ) 

P{T  +  C) 


3/4 


,  0-5148  0.2683 

1  -0.0618«-»- - -I- 


u 


(75) 


where  C  =  120  K,  T2  is  a  2-m  a.s.l.  temperature  measurement  (in  degrees  Kelvin),  and 
U2  is  a  2-m  a.s.l.  windspced  measurement  (in  m/s). 


3.  SUBROUTINE  DESCRIPTIONS 

The  TGRAD  program  is  designed  as  a  set  of  interactive  subroutines  that  are  based  around 
a  main  subroutine  C2illed  TGRAD.  The  particular  form  for  the  main  routine  is  of  little 
consequence  to  the  overall  operations  of  the  subroutines  and  will  not  be  considered  in 
depth  here.  The  reasoning  is  that  any  suitable  routine  that  passes  the  correct  information 
to  the  subroutines  may  be  used.  K.  Kimkel  (as  part  of  the  work  les^iing  to  Kunkcl  (1985)) 
has  produced  a  relatively  simple  driver  routine  that  is  listed  in  section  4. 


3.1  TGRAD 

TGRAD  is  the  main  subroutine  of  the  code.  It  is  the  only  routine  the  main  code  needs  to 
call.  The  name  TGRAD  refers  to  the  routine’s  estimation  of  temperature  gradient  that  was 
the  principal  purpose  of  the  original  code.  Since  the  majority  of  the  code  is  based  on 
Deardorff  (1978),  equations  based  on  Dcardorff  are  annotated  by  the  capital  D  specifier  in 
the  source  code  listing. 

The  input  to  TGRAD  usually  involves  data  entered  through  two  separate  files.  In  the 
first,  a  description  of  the  site  is  given.  This  description  should  not  vary  significantly 
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with  time.  It  includes  data  on  the  site’s  latitude,  the  soil  type  (such  as  loam,  sand, 
and  clay),  the  fractional  foliage  cover,  the  predominant  cloud  type,  the  height  at  which 
the  atmospheric  characteristics  (for  example,  temperature  gradient  and  C^)  should  be 
estimated,  the  surface  roughness  length,  and  the  degrees  of  longitude  the  site  is  west  of  the 
standard  meridian.  The  second  file  contains  the  meteorological  data  from  a  measurement 
station  (assumed  to  be  2  m  a.s.l.)  coupled  with  the  corresponding  time  information  (Julian 
date  and  local  time  for  each  weather  record).  The  required  input  data  are  the  windspeed 
in  m/s,  temperature  in  degrees  Celsius,  fractional  cloudcover,  pressure  in  millibars,  and 
relative  humidity  in  percent.  This  information  is  passed  to  TGRAD  routine  by  the  driver. 

When  the  routine  is  called  for  the  first  time,  a  flag  will  be  set  to  zero,  indicating  that  the 
routine  must  set  up  default  initializations  of  the  surface  energy  budget  peirameters.  If  it  is 
not  the  first  time  the  routine  is  called,  the  routine  will  compute  the  energy  budget  based 
on  the  new  inputs  and  the  variable  values  stored  in  the  SAVED  common  block  of  passed 
data.  (See  the  list  of  parameters  and/or  the  program  comments  for  the  descriptions  of 
the  variables  in  the  various  common  blocks.)  Normally  the  initialization  flag  will  be  set 
to  one,  indicating  the  new  data  is  consistent  with  what  has  been  provided  previously.  If 
gaps  occur  in  the  data,  a  zero  can  again  be  passed  to  the  TGRAD  routine  to  allow  a  new 
initialization.* 

Once  the  subroutine  TGRAD  has  set  up  initial  conditions  (if  necessary),  soil  parameters 
are  also  set.  Then,  the  times  of  sunrise  and  sunset  are  found  from  subroutine  SUN,  the 
vapor  pressure  and  specific  humidity  of  the  air  are  found  from  subroutine  SPEHU,  groups  of 
constants  are  found  that  are  used  later  in  the  program,  and  the  amount  of  incident  solar 
radiation  is  found  from  subroutine  CLOUD.  The  specifics  of  these  subroutines  are  discussed 
in  subsequent  sections. 

Table  1  provides  a  listing  of  the  variables  used  in  the  subroutine  and  a  brief  description 
of  each.  Table  2  provides  a  listing  of  the  inputs  required  to  the  subroutine  and  their 
description. 


*Such  gaps  have  occurred  when  the  model  was  used  to  predict  the  behavior  of  the 
surface  layer  atmosphere  from  inputs  of  weather  observation  data  provided  by  the  .\ir 
Force’s  ETAC  (Environmental  Technical  Applications  Center)  database  of  airport  weather 
observations.  Certeiin  sites  do  not  make  weather  observations  at  night  if  no  aircraft  are 
scheduled  to  leind.  Six  hours  is  about  the  longest  gap  permissable  without  incurring 
significant  estimation  errors.  As  should  be  obvious,  the  TGRAD  model  is  not  predictive  in 
the  temporal  sense.  That  is,  it  cannot  operate  without  a  continuous  input  of  observed 
station  data.  Rather,  it  is  designed  to  fill  in  knowledge  gaps  concerning  the  vertical 
structure  of  the  surface  layer  using  the  meeisurements  from  a  single  tower  level  and  some 
indication  of  initial  soil  characteristics.  However,  if  a  synoptic  estimation  of  the  predicted 
weather  were  avaulable  (temperature,  winds,  relative  humidity,  and  cloud  cover  as  functions 
of  time),  estimates  of  surface  layer  flux  amd  vertical  structures  could  be  predicted,  based 
on  the  outcomes  of  the  higher  level  model. 
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TABLE  1.  SYMBOL  TABLE  FOR  SUBROUTINE  T6RAD 


Name 

Type 

Location 

Description 

AA,  AAA 

Real 

Local 

Coefficients  of  sensible  heat 

ABl,  AB2 

Real 

Local 

Coefficients  used  in  foliage  temperature  calculations 

ALBEDO 

Real 

Local 

Original  soil  albedo 

ALPHAF 

Real 

Local 

Foliage  albedo 

ALPHAG 

Real 

Local 

Soil  albedo  with  variable  moisture  content 

ALPHAP 

Real 

Local 

Fractional  surface  moisture 

AHUM 

Double 

Local 

Numerator  of  ground  temperature  equation 

BB,  BBB 

Real 

Local 

Coefficients  of  latent  heat 

BETA 

Real 

Local 

Constant  in  turbulence  equations 

BG 

Real 

Local 

FVactional  amount  of  bare  ground 

Cl,  C2 

Real 

Local 

Constants  used  in  force-restore  rate  equation  for  ground 
surface  temperature 

CCl,  CC2 

Real 

Local 

Coefficients  in  the  rate  equation  for  ground  surface 
moisture 

CCV(3) 

Real 

Common 

Cloud  fractions  for  three  cloud  layers:  1-highest,  2- 
middle,  S-lowest 

CDEGK 

Double 

Local 

Conversion  from  Celsius  to  Kelvin 

COSZ 

Real 

Local 

Cosine  of  the  solar  zenith  angle 

CP 

Real 

Local 

Specific  heat  of  fur  at  constant  pressure 

D1 

Real 

Local 

Soil  depth  influenced  by  the  diurnal  temperature  cycle, 
equal  to  SqRT(THERDF*TAUi) 

D12,  DIG 

Real 

Local 

Values  of  Di  with  moisture  added 

DIP 

Real 

Local 

Soil  depth  (0.1m)  influenced  by  the  soil  moisture  cycle 

D2 

Real 

Local 

Soil  depth  influenced  by  the  annual  temperature  cycle, 
equal  to  19.1  xDl 

Name 

Type 

Location 

Description 

D2P 

Real 

Local 

Soil  depth  (0.5m)  influenced  by  the  seasonal  moisture 
variations 

DBARD2 

Real 

Local 

Earth-sim  distance  modification  factor 

DEC 

Real 

Local 

Declination  of  the  sun  in  radians 

DELC 

Real 

Local 

Step  function,  equal  to  1,  except  during  condensation 

DELLHG 

Real 

Common 

Difference  between  the  local  longitude  and  the  time  zone 
meridian  (west  >  0) 

DELTXK 

Real 

Local 

Diffence  between  TEFF  and  TG  x  K 

DEVON 

Double 

Local 

Denominator  of  ground  temperature  equation 

DLAT 

Real 

Common 

Local  latitude  in  degrees  (30.5  is  30°  30’) 

DQSGDT 

Real 

Local 

Derivative  of  QSG  with  respect  to  TG 

DT 

Real 

Local 

Time  increment  in  seconds 

DTDZ 

Real 

Common 

Temperature  gradient 

EF 

Real 

Local 

Evaporation  rate  from  foliage 

EFFALB 

Real 

Local 

Ground-foliage  effective  albedo 

EFPOT 

Real 

Local 

Potential  evaporation  rate  from  foliage 

EG 

Real 

Local 

Evaporation  rate  at  ground  surface 

EH 

Real 

Local 

Evaporation  rate  just  above  foliage  canopy 

Real 

Local 

Foliage  emissivity 

ENG 

Real 

Local 

Ground  surface  emissivity 

ETR 

Real 

Local 

IVanspiration  rate 

FC 

Real 

Loccd 

FVactional  foliage  cover 

G1,2.3,4 

Real 

Local 

Sets  of  values  used  more  than  once 

HA 

Real 

Local 

Sum  of  fluxes  to  atmosphere  (positive  when  directed 
upward) 

HETCAP 

Real 

Local 

Heat  capacity  of  soil 

HRAIGL 

Real 

Local 

Hour  angle  of  sim  in  radians 
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Name 


HSF 


HSG 


HSH 


ICL 


IGAP 


ISOIL 


ISOILO 


ISUMRS 


ISUMST 


ITG 


lYR 


JD 


JDO 


JSOIL 


LAM2 


LAMG 


LANGO 


LF 


HIM 


NIHDEL 


OBUKLV 


0ME3RD 


Type 

Location 

Description 

Real 

Local 

Sensible  heat  flux  from  foliage  (positive  when  directed 
upward) 

Real 

Local 

Sensible  heat  flux  from  ground  surface 

Real 

Conunon 

Sensible  heat  just  above  foliage  canopy 

Integer 

Common 

Cloud  types  coding  flag 

Integer 

Common 

Flag  to  reinitialize  after  gaps  in  data 

Integer 

Local 

Soil  type  local  to  TCHAD  (can  be  changed) 

Integer 

Local 

Previous  soil  type 

Integer 

Common 

Time  of  local  sunrise  in  minutes  from  midnight 

Integer 

Common 

Time  of  local  sunset  in  minutes  from  midnight 

Integer 

Common 

Flag  for  using  measured  ground  temperature 

Integer 

Common 

Last  two  digits  of  the  year 

Integer 

Common 

Julian  date 

Integer 

Local 

Previous  Julian  date 

Integer 

Common 

Initial  soil  type 

Real 

Local 

Von  Karman’s  constant 

Real 

Local 

Latent  heat  of  vaporization 

Real 

Local 

Thermal  conductivity  of  deep  soil 

Real 

Local 

Thermal  conductivity  of  top  soil 

Real 

Local 

Initial  thermal  conductivity  of  soil 

Real 

Local 

Latent  heat  of  fusion 

Integer 

Common 

The  current  minute  of  the  day  from  midnight 

Integer 

Common 

Time  increment  in  minutes 

Real 

Local 

Net  leaf  area  index 

Real 

Common 

Monin-Obukhov  length 

Double 

Local 

Constant  equal  to  one- third 
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Name 

Type 

Location 

Description 

PI 

Real 

Local 

Value  multiplied  by  HA  in  force  restore  equation  for 
grotmd  temperature  (D-8a) 

P2 

Real 

Local 

Value  multiplied  by  (TG-T2)  in  same  equation 

PC2 

Real 

Local 

Heat  capacity  of  deep  soil 

PCG 

Real 

Local 

Heat  capacity  of  top  soil 

PG 

Real 

Local 

Precipitation  rate  felt  at  groimd  siirface 

PI 

Reeil 

Local 

Constant  equal  to  tt 

PRECIP 

Real 

Common 

Precipitation  rate 

PRES 

Real 

Common 

Atmospheric  pressure  in  millibars 

PSCSDl 

Real 

Local 

Heat  capacity  times  thermal  diffusivity  for  top  soil 

PSCSD2 

Real 

Local 

Heat  capacity  times  thermal  diffusivity  for  deep  soil 

QA 

Real 

Local 

Specific  humidity  of  the  air 

QDUM 

Real 

Local 

Dummy  variable  used  in  subroutine  call 

qsA 

Real 

Local 

Saturated  specific  humidity  of  the  air 

qsF 

Real 

Local 

Saturated  specific  humidity  of  the  foliage 

qsG 

Real 

Local 

Saturated  specific  humidity  of  the  ground 

RA 

Real 

Local 

Aerodynamic  resistance  (returned  from  CN) 

RAIR 

Real 

Local 

Gas  constant  for  air 

RH 

Real 

Common 

Relative  humidity  (not  in  percent) 

RHOA 

Real 

Local 

Density  of  the  air 

RHOW 

Real 

Local 

Density  of  water 

RL 

Real 

Local 

Longwave  radiative  flux  from  the  atmosphere 

RLAT 

Real 

Local 

Local  latitude  in  radians 

RLG1,2 

Real 

Local 

Sets  of  values  of  RL  used  more  than  once 

RP 

Real 

Local 

Soil  moisture  interpolation  factor 

RPP 

Real 

Local 

FVaction  of  potential  evaporation  rate  from  foliage 
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Name 

Type 

Location 

Description 

RS 

Real 

Local 

Generalized  stomatal  resistance 

SO 

Real 

Local 

Maximtun  incident  solar  radiation 

SDI 

Real 

Common 

Incident  solar  radiation  at  the  ground  surface 

SGHEAT 

Real 

Local 

Net  solar  radiative  flux 

SIG 

Real 

Local 

Stefan-Boltzman  constant 

SNAX 

Real 

Local 

Maximiun  solar  radiation  at  locsJ  noon 

S0IL(9,6) 

Real 

Local 

Sets  of  soil  parameters 

SOLNAX 

Real 

Local 

Solar  constant 

SS 

Real 

Local 

Seasonal  dependence  of  stomatal  resistance 

T2 

Double 

Local 

Temperature  of  the  soil  over  depth  D2  in  Kelvin 

T2C 

Real 

Common 

Same  as  T2  in  **  Celsius 

TAG 

Real 

Common 

Air  temp^ature  at  reference  height  in  °  Celsius 

TAIR 

Double 

Local 

Air  temperature  at  reference  height  in  Kelvin 

TAUl 

Real 

Local 

Diurnal  period  in  seconds 

TCCV 

Real 

Common 

Total  fraction  of  cloud  cover 

TEFF 

Double 

Local 

Groimd  surface/foliage  effective  temperature 

TF 

Double 

Local 

Average  foliage  temperature  in  Kelvin 

TFC 

Real 

Common 

Average  foliage  temperature  in  “  Celsius 

TG 

Double 

Local 

Ground  surface  temperature  in  Kelvin 

TGC 

Real 

Common 

Ground  surface  temperature  in  **  Celsius 

THERDF 

Real 

Local 

Thermal  diffuavity  of  soil 

THOOI 

Real 

Local 

Time  of  local  noon  in  minutes  from  midnight 

TSTAR 

Real 

Local 

Virtual  temperature  in  Kelvin 

TUA 

Double 

Local 

Upper  air  temperature  in  Kelvin 

USTAR 

Real 

Common 

FViction  velocity 

VPA 

Real 

Local 

Vapor  pressure  of  the  air  in  millibars 

45 


Name 

Type 

Location 

Description 

VPDUM 

Real 

Local 

Dtimmy  variable  iised  in  subroutine  call 

W2 

Double 

Local 

Dimensionless  volumetric  concentration  of  soil  moisture 
over  depth  D2 

V20 

Double 

Local 

Initial  value  of  U2 

vn> 

Double 

Local 

Mass  of  liquid  water  retained  by  foliage  per  unit 
horizontal  groimd  area 

WDMAX 

Double 

Local 

Maximum  value  of  UD  beyond  which  runoff  to  the 
ground  occurs 

WG 

Double 

Local 

Same  as  U2  for  top  soil 

UGO 

Double 

Local 

Initial  value  of  UG 

WIND 

Real 

Common 

Windspeed  at  reference  height 

UK 

Double 

Local 

Critical  or  saturated  value  of  UG 

UMAX 

Real 

Local 

Maximum  value  of  UG 

WMIU 

Real 

Local 

Minimum  value  of  UG 

US 

Real 

Local 

0.9U2  +  O.lUG 

UUILT 

Real 

Local 

Wilting  point  value  of  UG 

Z 

Real 

Common 

Desired  height  for  temperature  gradient 

ZO 

Real 

Common 

Roughness  height  of  the  terrain 

TABLE  2.  INPUTS  TO  SUBROUTINE  TGRAD 


Name 

Description 

^-indicates  the  value  has  been  preset  within  the  subroutine 

lYR 

Year  used  in  calculating  solar  declination. 

JD 

Julian  date. 

JS0IL,S0IL« 

JSOIL  indicates  soil  type,  SOIL  is  an  array  containing  values  of  soil 
characteristics  arranged  by  soil  type.  The  following  are  elements  of 
SOIL. 

SOILd,*)* 

Ground  surface  emissivity.  Set  to  EMG  in  main  routine. 

S0IL(2,*)* 

ALBEDO  of  the  surface. 

S0IL(3,*)* 

THERDF  (thermal  difhisivity)  of  the  ground. 

S0IL(4,*)* 

HETCAP  (heat  capacity)  of  the  ground. 

S0IL(5,*)* 

HGO,  the  initial  fractional  water  content,  by  volume,  of  the  first  10  cm 
of  soil. 

S0IL(6,*)* 

WK,  the  fractional  wat^  content,  by  volmne,  at  which  the  soil  acts  as 
if  satxirated. 

S0IL(7,*)* 

W20,  the  initial  fractional  water  content,  by  volume,  of  the  first  50  cm 
of  soil.  The  value  of  VGO  is  subsumed  within  this  value. 

S0IL(8,*)* 

WWILT,  is  a  measure  of  the  resistance  of  foliage  to  evaporative  processes 
as  the  soil  becomes  drier. 

S0IL(9,*)* 

VNII,  the  minimum  value  allowed  for  the  soil  moisture.  If  no  minimtun 
was  given,  the  calculation  of  VS  could  reach  zero  and  the  stomatal 
resistance  equation  would  approach  infinity. 

MIX 

The  driver  routine  accoimts  for  time  and  reads  data  from  an  external 
sovirce  on  an  hourly  to  3  hourly  basis.  Mil  is  the  minutes  past 
midnight  in  local  standard  time. 

MIIDEL 

Time  increment  between  successive  surface  energy  budget  calcula¬ 
tions.  A  standard  value  of  10  minutes  has  been  chosen. 

ITG 

ITG  flags  whether  or  not  user  wants  to  input  soil  temperature. 
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Name 

Description 

IGAP 

IGAP  is  a  driver  flag  indicating  that  the  time  since  the  last  data  set 
was  too  long  to  be  relevant.  Parameters  are  reinitialized  as  at  the 
start  of  the  program. 

ISOILO 

ISOILO  is  the  initial  soil  type;  provision  is  made  within  TGRAD  to 
vary  the  soil  type  by  season,  so  in  cases  of  reinitialization,  the  original 
soil  parameters  can  be  recalled  if  needed. 

ICL 

Cloud  type  category,  as  listed  in  the  CLOUD  subroutine. 

PRES 

Measured  pressure,  in  millibars. 

RH 

Measured  relative  humidity. 

WIHD 

Measured  windspeeds  converted  to  meters  per  second. 

CCV(3) 

Cloud  cover  amounts  in  three  cloud  height  categories:  low,  medium, 
and  high.  Only  occasionally  are  three  layers  of  data  available. 

TCCV 

Total  fractional  cloud  cover. 

DLAT 

Latitude  of  site  chosen. 

AFC 

The  original  fractional  foliage  cover  value,  used  when  reinitializing. 

ZO 

Surface  roughness  length. 

DELLNG 

DELLNG  is  the  difference  in  longitude  between  the  site’s  longitude  and 
the  meridian  used  as  a  basis  of  time.  For  example,  WSMR  is  located 
at  approximately  106®  W  longitude.  It’s  DELLNG  would  be  +1®.  A 
location  in  the  eastern  hemisphere,  at  106®  E  would  have  a  DELLNG 
of  -1®.  This  convention  is  used  to  calculate  the  exact  time  of  sunrise 

and  sunset  for  each  location. 

Z 

Height  desired  for  temperature  gradient  calculation.  According  to 
similarity  theory,  knowledge  of  certain  scaling  parameters  allows  one 
to  calculate  the  gradient  anywhere  within  the  boundary  layer. 

PRECIP 

Precipitation  rate  in  kg/s-m*. 

TAG 

Measured  station  temperature  in  ®  Celsius. 

TGC 

Ground  temperature  in  ®  Celsius.  While  usually  calculated  within 
the  model,  it  can  also  be  entered  as  a  program  input  variable. 
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Name 

Description 

K* 

Von  Karman’s  constant,  a  similarity  theory  constant  equal  to  0.4. 

ALPHAF« 

Albedo  of  the  foliage,  set  to  0.2. 

EMF* 

Emissivity  of  the  foliage.  Albedo  is  a  general  refiectivity  to  the 
solar  radiation  spectrum,  and  emissivity  is  the  generalized  ability  of 
a  surface  to  emit  energy  in  the  infrared  spectrum  through  ‘graybody’ 
radiation. 

SOLMAX* 

The  solar  constant  of  1369.2  W/m^  at  the  Earth’s  orbit  radius. 

ALBSN0« 

Albedo  of  a  typical  snow  cover,  set  at  0.6. 

BETA* 

A  similarity  constant  used  for  nocturnal  boundary  layer  character¬ 
ization.  In  the  literature  varies  between  approximately  4.7  and 
9.5. 

The  soil  parameters  are  declared  in  the  BLOCK  DATA  segment  of  the  code.  Deardorff 
(1978)  listed  5  sets  of  soil  parameters  in  his  paper.  Other  sets  were  found  in  Oke  (1978). 
These  latter  sets  had  separate  categories  for  dry  and  wet  soil  types,  while  Deardorff  had 
simple  parameterizations  of  albedo,  heat  capacity,  and  thermsd  conductivity  to  account 
for  variations  in  soil  moisture.  Deardorff’s  equations  were  adopted  as  equations  (54a) 

through  (54c),  except  that  (54b)  has  been  correctly  modeled  by  including  the  wUlx  factor 
in  the  denominator.  A  comparison  between  Deardorff’s  soils  and  Oke's  soils  shows  they 
are  virtually  the  same  soil  types  when  moisture  differences  are  taken  into  accoimt.  The 
values  of  Wk  and  were  taken  from  a  chart  foimd  on  page  42  of  Oke  (1978).  Soil  types 
include  clay,  snow,  a  “wet”  desert  sand  (steppe),  and  a  dry  desert  sand.  Alternatively, 
different  soil  parameters  could  be  read  in  for  the  location,  but  this  would  require  a  change 
in  the  parameters  passed  through  the  unnamed  common  block. 

3.2  CN 

Subroutine  CN  uses  the  theory  in  section  2.1.4  to  calculate  various  flux-profile  parameters, 
including  the  Obukhov  length,  OBUKLN;  the  friction  velocity,  USTAR;  and  the  dimensionless 
temperature  gradient  function,  PHIH.  PVom  these  quantities  the  sensible  heat  fiux  and 
the  other  propagation  related  parameters  discussed  in  section  2.4  can  be  calculated.  In 
the  program,  equations  derived  from  Thom’s  work  are  designated  with  a  T.  The  basic 
flux-profile  equations  are  those  obtained  from  Haima  ti  al.  (1982).  These  equations  are 
designated  with  an  H,  though  several  other  authors  have  identi^  equations.*  Table  3 
shows  the  names  and  short  descriptions  of  all  the  variables  used  in  subroutine  CN. 


*Also,  for  lack  of  characters  at  least  one  variable  needed  to  be  renamed.  r//,  the 
atmospheric  resistance  term,  is  expressed  as  RA,  while  Rh,  the  relative  hximidity  is 
expressed  as  RH. 
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TABLE  3.  SYMBOL  TABLE  FOR  SUBROUTINE  CN 


Name 

Type 

Location 

Description 

Bi 

Real 

Local 

Excess  of  resistance  due  to  heat  flux 

CHARGE 

Real 

Local 

Relative  change  in  USTAR  iterations 

CON 

Real 

Local 

A  group  of  constants 

DELZ 

Real 

Local 

Desired  height  (Z)  -  roughness  height  (ZO) 

DELZS 

Real 

Local 

Reference  height  (STAHT)  -  roughness  height 

DPTDZ 

Real 

Local 

Potential  temperature  gradient 

HS 

Real 

Local 

Sensible  heat  flux/(RH0A*CP) 

H 

Integer 

Local 

DO  loop  counter 

KH 

Real 

Local 

Eddy  diflusivity  for  heat 

OBMIM 

Real 

Local 

Minimiun  value  of  OBUKLN  used  at  night 

PHIH 

Real 

Local 

Dimensionless  temperature  gradient 

PHIHMX 

Real 

Local 

Maximtun  allowed  value  of  PHIH  at  night 

PHIN 

Real 

Local 

Dimensionless  wind  gradient 

POH 

Real 

Local 

Exponent  used  in  conversion  from  potential  tempera- 
ttire  to  temperature 

PSI 

Real 

Local 

The  result  of  integrating  PHIN 

STAHT 

Real 

Local 

Reference  height  for  air  temperatmre  and  windspeed 
measurements 

UST 

Real 

Local 

Previous  value  of  USTAR 

USTMII 

Real 

Local 

Minimum  value  of  USTAR  used  at  night 

3.3  CLOUD 


Subroutine  CLOUD  was  drawn  almost  verbatim  firom  Shapiro’s  (1982)  three  layer  cloud 
model.  This  is  an  empirical  model  based  on  statistical  averages  only.  Averages  represent 
mean  incident  radiation  from  the  sun  arriving  at  the  surface  of  the  earth.  These  mean 
results  average  out  the  fluctuations  that  would  occur  when  direct  sunlight  hits  the  surface 
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on  a  cloudy  day  versus  when  direct  sunlight  is  obscured  by  clouds.  Therefore  the  results 
of  the  model  equate  to  an  area  average  rather  than  a  point  calculation. 

The  parameters  brought  into  the  subroutine  include  the  cosine  of  the  zenith  angle  (COSZ), 
the  solar  constant  (SO),  and  the  effective  albedo  (equation  (15a))  of  the  foliage/groimd 
surface.  The  fractional  amoimts  of  clouds  in  the  three  layers  and  their  types  are  passed 
via  the  tmnamed  common.  The  only  output  from  the  routine  is  the  incident  shortwave 
flux  at  the  groimd.  This  result  does  not  distinguish  between  direct  and  diffuse  radiation. 

Table  4  is  a  listing  of  the  parameters  used  in  the  subroutine. 

The  variable  ICLQUD  used  in  the  subroutine  is  a  coded  number  used  to  denote  the  cloud 
types  in  each  of  the  three  layers.  This  is  accomplished  by  designating  individual  digits 
of  a  five  digit  number  to  zeros,  ones,  or  twos.  Each  digit  characterizes  one  of  the  three 
cloud  layers  or  designates  a  rain  or  fog  condition.  In  layer  one,  there  is  a  distinction 
made  between  thick  cirrus/cirrostratus  (Ci/Cs)  and  thin.  Layer  two  does  not  distinguish 
between  alto-type  clouds.  Layer  three  distinguishes  between  stratus  and  cumulus  tj'pes 
and  also  allows  for  fog  and/or  smoke.  Rain  is  denoted  by  overcast  conditions  in  all  layers 
and  thick  Ci/Cs  in  layer  one.  The  fractional  amounts  of  cloud  type  by  layer  are  given  in 
CCV(3)  and  are  reassigned  to  FI,  F2,  and  F3  within  the  subroutine.  Normally  only  one 
cloud  type  is  set  in  this  integer  because  the  input  data  available  usually  does  not  specify 
cloud  type.  However,  inclusion  of  a  cloud  type  is  a  simple  addition. 

Except  for  the  following  notational  departures,  Shapiro’s  (1982)  work  was  faithfully 
reproduced.  The  overcast  coefficients  are  denoted  by  P  (instead  of  RHO)  and  U  (instead  of 
TAU),  and  one  character  is  used  for  Shapiro’s  lower-case  variable  names,  and  two  characters 
are  used  for  Shapiro’s  upper-case  variable  names.  The  reader  interested  in  more  detail 
regarding  Shapiro’s  model  is  turged  to  obtain  his  report.  His  methodology  and  findings 
will  not  be  repeated  here. 
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TABLE  4.  SYMBOL  TABLE  FOR  SUBROUTINE  CLOUD 


Name 

Type 

Location 

Description 

BH 

Real 

Local 

Excess  of  resistance  due  to  heat  flux 

Dl,2,3 

Real 

Local 

dl,  d2,  d3  from  Shapiro 

DD1,2 

Real 

Local 

D1  and  D2  from  Shapiro 

FI. 2, 3 

Real 

Local 

Same  as  CCV(3)  above 

ICLOUD 

Integer 

Local 

Same  as  ICL  above 

ID2,3 

Integer 

Local 

Flags  for  using  diffuse  values  for  layers  2  and  3 

IPU1,3 

Integer 

Local 

Flags  for  type  of  cloud  to  use  in  layers  1  and  3 

IRAIN 

Integer 

Local 

Flag  for  rain 

IRT3 

Integer 

Local 

Flag  for  fog  and/or  smoke  present  in  layer  3 

PI, 2, 3 

Real 

Local 

Overcast  reflectivities  for  3  layers 

PHI1,2,3 

Real 

Local 

Weighting  frmctions  for  3  layers 

Rl,2,3 

Real 

Local 

Clear  air  reflectivities  for  3  layers 

RG 

Real 

Formal 

Ground  stirface-foliage  albedo 

RR1,2,3 

Real 

Local 

Total  reflectivities  for  3  layers 

Tl,2,3 

Real 

Local 

Clear  air  transmissivities  for  3  layers 

TT1,2,3 

Real 

Local 

Total  transmissivities  for  3  layers 

Ul,2,3 

Real 

Local 

Overcast  transmissivities  for  3  layers 

Wl.2.3 

Real 

Loc2d 

Weighting  frmctions  for  3  layers 

XO 

Real 

Local 

Incident  solar  radiation  above  clouds 

X3 

Real 

Formal 

Incident  solar  radiation  at  groimd  siuface 
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3.4  TFOL 


The  computation  of  the  foliage  temperature  is  accomplished  in  subroutine  TFOL  using 
the  iterative  Newton-Raphson  procedure  described  in  section  2.2.  Newton’s  method  is  the 
fastest  method  for  solving  a  polynomial,  as  long  as  the  first  estimate  is  reasonably  close  to 
the  solution.  This  is  virtually  assured  by  the  small  change  expected  in  foliage  temperature 
between  timesteps.  Table  5  itemizes  the  variables  used  in  this  subroutine. 

TABLE  5.  SYMBOL  TABLE  FOR  SUBROUTINE  TFOL 


Name 

Type 

Location 

Description 

CHANGE 

Double 

Local 

Relative  change  of  TF  after  iteration 

DQSFDT 

Real 

Local 

Derivative  of  saturated  specific  humidity  of  foliage  with 
respect  to  TF 

FOFTF 

Double 

Local 

F(TF) 

FPOFTF 

Double 

Local 

F»(tf) 

I 

Integer 

Local 

DO  loop  cotmter 

QDUM 

Real 

Local 

Dummy  variable  used  in  subroutine  call 

VPDUM 

Real 

Local 

Dummy  variable  used  in  subroutine  call 

In  each  iterative  step,  the  saturated  specific  humidity  of  the  foliage  and  its  derivative  are 
found  by  calling  subroutine  SPEHU.  Then  FOPTF  (meaning  the  function  F  o{0j)  and  FPOFTF 
(F'  in  equation  (47))  are  found  and  the  change  in  TF  {Bj)  is  found.  TF  is  changed  by  this 
amount  and  a  check  is  made  on  the  amount  of  the  change.  A  decision  was  made  to  exit 
the  subroutine  when  the  change  is  less  than  10“^  K.  If  convergence  does  not  occur  within 
20  iterations,  the  program  writes  an  error  message  along  with  the  last  value  of  TF.  The 
program  has  never  encountered  a  problem  attempting  to  estimate  Bf. 

3.5  SUN 

Subroutine  SUI  finds  the  times  of  local  sunrise,  simset,  and  noon  (sun  directly  south)  as 
well  as  the  solar  declination.  It  requires  the  Julian  date  ( JD),  the  latitude  in  radians  (RLAT), 
and  the  difference  in  degrees  between  the  location  and  the  standard  meridian  used  for  the 
local  time  (DELLHG).  The  latitude  enters  TGRAD  through  the  unnamed  common  as  DLAT 
(the  latitude  in  degrees).  Then,  RLAT  is  calculated  and  passed  from  TGRAD  to  SUN.  Both 
DLAT  and  DELLHG  are  in  decimal  degrees.  The  subroutine  is  used  to  compute  equations  (7) 
through  (12),  which  are  used  to  set  up  equation  (14),  the  equation  for  the  zenith  angle. 
Table  6  is  a  compendium  of  the  different  variables  used  in  subroutine  SUI. 
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TABLE  6.  SYMBOL  TABLE  FOR  SUBROUTINE  SUN 


Name  Type  Location  Description 

AHOOH  Real  Local  Solar  angle  at  local  noon 

BETA  Real  Local  Definition  of  angle  of  stmrise  or  sunset 

EOT  Real  Local  Value  from  equation  of  time 

HRAMGL  Real  Local  Hour  angle  of  the  sim 

SBETA  Real  Local  SIN  (BETA) 


3.6  SPEHU 

Subroutine  SPEHU  takes  a  temperature,  a  relative  humidity,  and  the  local  pressme  and 
calculates  vapor  pressure,  specific  hiunidity,  saturated  specific  humidity,  and  the  derivative 
of  the  saturated  specific  humidity.  There  are  a  total  of  three  calls  to  this  routine;  two 
from  TGRAD  and  one  from  TFOL. 

Equations  from  Haurwitz  (1941)  were  used,  but  these  equations  can  be  found  in  most 
meteorological  texts.  The  derivative  of  the  saturated  specific  humidity  with  respect  to 
temperature  is  a  straight  forward  differentiation,  performed  in  appendix  A. 

TABLE  7.  SYMBOL  TABLE  FOR  SUBROUTINE  SPEHU 


Name 

Type 

Location 

Description 

DQSDT 

Real 

Formal 

Derivative  of  the  saturated  specific  humidity  with 
respect  to  the  temperatme 

POW 

Real 

Local 

Exponent  used  in  vapor  pressure  equation 

q 

Real 

Formed 

Specific  humidity 

qs 

Real 

Formal 

Satmated  specific  humidity 

TC 

Real 

Formal 

Temperatme  in  ®  Celsius 

VP 

Real 

Formal 

Vapor  pressxue 

VPS 

Real 

Local 

Satmated  vapor  pressure 
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4.  SAMPLE  DRIVER 


As  noted,  subroutine  TCHAD  is  called  by  a  driver.  It  was  written  this  way  to  accommodate 
differences  in  the  formats  of  various  meteorological  data  bases.  Originally,  the  Ktmkel  code 
was  modified  to  deal  with  two  different  sources  of  data,  but  even  further  applications  were 
foimd  to  be  possible  in  the  modular  form  adopted.  The  sole  changes  required  were  in  the 
driver.  Generally  this  is  true,  and  since  a  user  will  likely  have  a  unique  application  and 
data  source,  the  driver  listed  as  figure  3  will  serve  as  a  model  for  a  variety  of  applications. 

Figure  4  shows  a  sample  data  file.  The  data  was  made  available  from  meteorological 
observations  taken  at  C  Station,  White  Sands  Missile  Range,  New  Mexico.  Windspeeds 
recorded  as  zeros  were  adjusted  upwards  by  0.5  m/s  to  avoid  computational  errors  within 
the  CN  subroutine.  The  Julian  date  refers  to  2  November  1984.  When  running  the 
program,  it  will  be  necessary  to  remove  the  first  line  from  the  following  table.  The  header 
line  is  provided  as  a  simplifying  reference  to  the  necessary  inputs. 

5.  MODEL  PERFORMANCE  AND  SENSITIVITY  ANALYSIS 

The  performance  of  the  TGRAD  model  has  been  evaluated  twice.  First,  Kunkel  (1985) 
evaluated  the  ability  of  the  model  to  estimate  the  vertical  temperature  gradient.  Then 
Tofsted  and  Gillespie  (1986)  compared  the  model  performance  with  the  Deardorff  model 
and  an  improved  version  of  TGRAD,  which  we  shall  call  TGRAD  2.0.  Kunkel’s  analysis 
was  the  motivation  for  TGRAD  2.0,  since  his  analysis  uncovered  two  areas  where  TGRAD 
performed  below  expectations.  The  first  area  was  in  regard  to  using  a  surface-energy- 
budget-driven  near-surface  air  temperature.  As  indicated  in  section  2.1,  this  modification 
has  been  included  in  this  version  of  the  model.  A  second  area  was  model  performance 
for  snow-covered  terrain.  A  follow-on  version  of  TGRAD,  called  TGGEN,  attempted  to 
account  for  this  effect  through  inclusion  of  melting  effects.  However,  a  decision  was  made 
to  document  the  TGRAD  model  without  describing  the  snow  interaction  effects  since  the 
major  portion  of  the  model  has  not  changed  and  because  the  amount  of  description  required 
to  explain  the  snow-cover  effects  is  rather  large  and  best  left  to  a  separate  discussion. 

Kunkel's  other  basic  findings  related  to  the  direct  components  in  the  equation  computing 
temperature  gradient.  These  factors  were;  the  degree  of  accuracy  in  the  measurement 
of  windspeed  and  temperature,  the  model's  method  used  in  estimating  sensible  heat, 
error  in  estimating  the  roughness  length  at  the  site,  and  differences  in  the  theoretical 
functional  form  for  the  <f>|^  term.  He  found  first  that  the  measurement  of  temperature  was 
unimportant.  Second,  he  determined  that  changes  in  the  roughness  length  of  50  percent 
only  produced  gradient  changes  up  to  16  percent,  indicating  relative  model  insensitivity 
to  large  errors  in  estimating  the  roughness  length.  Third,  he  analyzed  the  ability  of  the 
model  to  estimate  the  sensible  heat  flux.  He  foimd  the  model  does  a  reasonable  job  for  a 
variety  of  different  surface  types.  His  table  was  left  out  of  the  proceedings  of  his  EOSAEL 
conference  paper,  so  it  is  included  here  as  table  8.'*' 


*  As  Kunkel  notes,  the  values  available  in  the  literature  are  few,  and,  in  general,  insufficient 
information  is  provided  with  the  sensible  heat  flux  values  to  be  able  to  accurately  use  the 
model.  For  this  reason  Kunkel  ran  the  model  with  a  range  of  inputs  in  order  to  bracket 
the  possible  conditions  present  when  the  cited  flux  values  occurred. 
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PROGRAM  TGS  ! 1000  <860203. 0926> 

C 

C  DRIVER  PROGRAM  FOR  TGRAD  SURFACE  EHERGV  BUDGET  MODEL 
C 

INTEGER  PP(5) 

REAL*8  MG, UMAX 

CraOION  JVR,  JD,  JSOIL, JMIM ,  JMIMDEL,  JTG ,  J6AP ,  JSUMRS ,  JSUNST ,  ISOI LO , 

1  JCL,PRES,RH,MIND,CCV(3),TCCV,DLAT,DELLNG,AFC,20,Z,PRECIP,DTDZ, 

2  SDN,HSH,TAC,TGC,TFC.T2C,OSTAR,OBUiaM.HG,MMAX,CN2 
C 

C  FF(1)  IS  THE  ENDING  MINUTE  FROM  MIDNIGHT  OF  THE  JULIAN  DATE 

C  FROM  UNIT  4t  CAN  BE  GREATER  THAN  1440. 

C  PF(3)  IS  THE  OUTPUT  INTERVAL  FOR  UNITS  9  AND  10. 

C  PP(3)  IS  THE  STARTING  MINUTE  FROM  MIDNIGHT  OF  THE  JULIAN  DATE 

C  FROM  UNIT  4. 

C  PP(4)  IS  EITHER  THE  VALUE  OF  THE  TIME  INTERVAL  USED  IN  TGRAD, 

C  (ALLOWABLE  VALUES:  1-10),  OR  THE  VALUE  OF  THE  RANGE  IN 

C  THE  MISS  DISTANCE  CALCULATIONS  (ALLOWABLE  VALUES;  >500  M)  . 

C  PP(S)  IS  A  FLAG  FOR  USING  THE  MEASURED  GROUND  TEMPS; 

C  DEFAULT  IS  OFF. 

C 

CALL  RMPAR(PP) 

C  1440  30  0  10  0 
INT-PP(2) 

IF(PP(2) .LE.0)INT>1 
IENC>>PP(1) 

IF ( lEND. LE . 1) IEND>32767 
ISTART-PP(3) 

IF(ISTART.  LT.  0)  ISTART-1440+ISTAHT 
ISTART>MOD(ISTART.  1440) 

ADT>0.0 

JMIH0EL*10 

RANGE-1000.0 

IF(PP(4) .GE.S00)RANGE-PP(4) 

JTG-PP(5) 

JSTAT-0 

HRITE(9, 700) 'TGRAD  RESULTS' 

700  FORMAT (///////A13/) 

REWIND  4 
REWIND  8 
JGAP-0 
IBA1>0 
JHIN-ISTART 
C 

C  READ  IN  LOCATION  PARAMETERS;  MORE  THAN  ONE  SET  CAN  BE  USED  FOR 
C  THE  SAME  MET  DATA. 

C 

READ ( 4 , • ) JDO , AFC , Z , Z  0 , DLAT , DELLMG , ICLOUD , JSOIL , JYR , PRETRG 
JD-JDO 

HRITE(7,*)JD, AFC, Z,Z0, DLAT, DELUIG, ICLOUD, JSOIL, JYR, IWHICH 
GO  TO  100 
C 
C 

20  IF(MOD(JMIN->1440,JMINDEL)  .NE.O)CO  TO  100 
I F  ( JGAP .  EQ .  0 .  OR .  JTG .  EQ .  2  )  TGOATG 
PRECIP-0.0 

IF (TAC. LE. PRETRG) PRECIP-2 . 77778E-4 
IF (JGAP. EO.O) THEN 
C 

C**  NEW  DAY:  BEGINNING  OF  DATA  OR  GAP  IN  DATA:  INITIALIZE  ** 

C 

HRITE(7,600) 

600  FORMAT  ('NEW  DAY') 

PI-3.14159 

RLAT-DUT*PI/1B0.0 

CALL  SUN(PI,RLAT,DEC,TNOON) 

ITIM1-JSUNR8/60*100'»NOD(JSUNRS,  60) 

ITIN2-JSUNST/60*1004MOD(JSUNST,  60) 

WRITE(9,») ' SUNRISE- ',ITIH1,'  SUNSET- ', ITIM2 
MRITE(9,690)  '  JO  MIN  WIND  TCCV  DTDZ  CN2  L  TG', 

1  '  TF  TA  SDN  N3/NNAX  IG  SNOW', 

2  '  HSH  MISSD  RADBAL' 

690  FORMAT(/A41,A38,A33/) 

JDl-JD 

ENDIF 

C 

C  NEW  JULIAN  DATE 

C 

IF(JD.NE.JDO)TKEN 

JDO-JD 

ITIMl-J8Ums/60»100*MOD(JSUHKS.  60) 

ITIM3-JSUNST/60*1004NOD(JSUNST, 60) 

MHITE(9,4) ' SUNRISE- ',ITIN1,'  SUNSET- ', ITIM2 
ENOIF 

Figure  3.  Simplified  driver  program  TGS. 
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c 

c 

CALL  TGItAO(JSTAT,WGAItF,AIS,ItAOB) 

AiaSSC-7 . 75E-S*Pi(ES/  (TAC-t'273 . 1C)  ••2*OTDZ*IIAllGB*RAMGE/2 . 0 
IP(lfOO(JMlM+1440,im)  .EQ.O)THBM 

NRITE ( 9 , 6 50 ) JD , JMIM , NIND . TCCV , DTDZ , CM2 , OBOKUl , 

1  TGC,TPC,TAC,SDM, 10000. •WGABP.HSH.AMISSC.RAOB 

650  FOIIIIAT(I4,I5,F6.2,F4.1,r5.2.r6.2,P7.0,3F5.1.P6.1,F9.4,12X, 

1  F9.4,F7.3,F7.2) 

ENDIF 

JiaM>JMlM+jmMOEL 
IF(JMIM.BQ.(IEMO»JlllMDEL))fiOTO  80 
IF(IFBItK(IO0M))80,70  11000 
70  IF(li0D(JinN,IMT).B().0)G0  TO  100 

IF((JD-JD1)*1440>JNIM.LT.IEMD}GO  TO  20 
100  REAO(a,*)ID,IT,PltE8,TAC,IIB,HIND,TCCV 

HRITE(7.*)JOO.JD,ID,IT,PllES,TAC,ItH,WIMD,TCCV 
IF ( ID. GT. JO) THEM 
JNIN-0 
JD-ID 
ENDIF 

IF(WIND.LT.0.5)NIND-0.5 
IF(in.  LT.  1. 0)  RH«RH*100.0 
ItH>RH/100.0 
JCLfOIIII 

IF  ( ICLOUO .  EQ .  2  )  JCL-02 11 1 

IF(ICL0UO.EQ.5)JCI^01121 

IF(ICL0UD.Ea.C)JCL-01112 

IF  ( ICtXXlD .  EQ .  7 )  JCI^12111 

CCV(l)-0.0 

CCV(2)«0.0 

CCV(3)-0.0 

I F ( ICLOUO . LE . 2 ) CCV ( 1 ) >TCCV 
IF ( ICLOOD . EQ . 3 ) CCV ( 2 ) >TCCV 
IF (ICLOUD. EQ . 4 .OR. ICLOUD.EQ. 5) CCV( 3 ) -TCCV 
GO  TO  20 
C 

C  AT  END  OF  FILE  OR  JD  IS  0 
C 

80  STOP 
END 

Figure  3.  (coni)  Simplified  driver  program  TGS. 


As  discussed  above  and  as  seen  in  table  8,  the  model  does  not  do  well  in  estimating  sensible 
heat  fiux  in  the  warm-advection-over-snow  case.  The  estimation  of  roughness  length  and 

differences  in  form  for  the  <f>h  function  were  the  most  significant  sources  of  error  for  this 
case. 

Kunkel  also  found  that  a  15  percent  error  in  measurement  of  the  windspeed  caused  at 
least  a  15  percent  change  in  the  estimated  temperature  gradient  at  night.  The  change 
WM  greatCTt  for  windspeeds  of  2  m/s  where  the  variation  was  24  percent.  Windspeeds  of 
this  velocity  are  common  for  highly  stable  episodes,  and  thus  windspeed  can  be  a  critical 
parameter  in  estimation  of  the  nocturnal  flux  cases. 

In  regard  to  Kunkcl’s  investigation  of  the  effects  of  variations  in  the  functional  forms 
for  <I>H,  Yaglom  (1977)  provided  a  compendium  of  various  functional  forms  proposed  by 
researchers  in  this  area.  Apparently,  even  after  twenty  years  of  research  in  the  area,  there  is 
no  general  agreement  on  the  proper  form  (especially  at  night).  In  his  analysis,  Kunkel  ran 
the  model  using  a  range  of  these  proposed  form-fit  coeflBcients.  The  results  showed  that 
for  the  higher  windspeeds  (between  10  and  15  m/s),  the  degree  of  variation  of  the  results 
about  the  standard  model  were  12  percent  during  the  day  and  17  percent  at  night.  At  6 
m/s  the  variations  were  20  percent  during  the  day  and  15  percent  at  night.  At  2  m/s  the 
variations  were  45  percent  during  the  day  and  57  percent  at  night.  Apparently  then,  the 
largest  source  of  error  is  due  to  the  functional  form  used  for  the  dimensionless  temperature 
gradient.  The  best  functional  form  and  related  coefficients  to  address  this  problem  are 
still  not  known.  Numerous  investigators  have  been  tackling  this  problem  for  the  past 
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JULD  TIME 

PRES 

TEMP 

RH 

WIMD 

CLCV 

226 

0 

880.0 

17.90 

73.0 

1.1 

0.00 

226 

30 

880.0 

18.18 

71.0 

1.1 

0.00 

226 

60 

880.0 

18.25 

70.0 

1.1 

0.00 

226 

90 

880.0 

17.10 

72.0 

1.1 

0.00 

226 

120 

880.0 

17.24 

76.0 

0.9 

0.00 

226 

150 

880.0 

17.43 

74.0 

0.9 

0.00 

226 

180 

880.0 

16.29 

76.0 

1.6 

0.00 

226 

210 

880.0 

15.47 

75.0 

1.5 

0.00 

226 

240 

880.0 

15.99 

73.0 

1.6 

0.00 

226 

270 

880.0 

16.79 

72.0 

1.3 

0.00 

226 

300 

880.0 

16.20 

71.0 

1.2 

0.00 

226 

330 

880.0 

15.44 

76.0 

1.3 

0.00 

226 

360 

880.0 

16.69 

74.0 

1.1 

0.00 

226 

390 

880.0 

17.96 

70.0 

1.4 

0.00 

226 

420 

880.0 

20.17 

59.0 

1.8 

0.00 

226 

450 

880.0 

21.71 

52.0 

1.6 

0.00 

226 

480 

880.0 

23.26 

46.0 

1.0 

0.00 

226 

510 

880.0 

24.73 

40.0 

0.8 

0.00 

226 

540 

880.0 

24.87 

39.0 

1.4 

0.00 

226 

570 

880.0 

25.88 

38.0 

2.1 

0.00 

226 

600 

880.0 

28.75 

33.0 

1.7 

0.00 

226 

630 

880.0 

30.32 

26.0 

1.5 

0.00 

226 

660 

880.0 

32.38 

21.0 

1.7 

0.00 

226 

690 

880.0 

33.25 

18.0 

2.3 

0.00 

226 

720 

880.0 

31.91 

15.0 

2.3 

0.00 

226 

750 

880.0 

34.27 

14.0 

2.8 

0.00 

226 

780 

880.0 

33.38 

12.0 

2.8 

0.00 

226 

810 

880.0 

33.01 

12.0 

4.5 

0.00 

226 

840 

880.0 

32.85 

11.0 

4.3 

0.00 

226 

870 

880.0 

33.14 

11.0 

3.8 

0.00 

226 

900 

880.0 

33.19 

11.0 

3.8 

0.00 

226 

930 

880.0 

33.11 

11.0 

3.6 

0.00 

226 

960 

880.0 

32.97 

11.0 

4.1 

0.00 

226 

990 

880.0 

32.70 

11.0 

3.9 

0.00 

226 

1020 

880.0 

32.55 

12.0 

3.9 

0.00 

226 

1050 

880.0 

31.71 

12.0 

4.3 

0.00 

226 

1080 

880.0 

30.54 

13.0 

3.2 

0.00 

226 

1110 

880.0 

28.90 

16.0 

2.3 

0.00 

226 

1140 

880.0 

26.66 

20.0 

1.6 

0.00 

226 

1170 

880.0 

25.13 

27.0 

1.4 

0.00 

226 

1200 

880.0 

24.30 

32.0 

1.6 

0.00 

226 

1230 

880.0 

22.71 

36.0 

1.7 

0.00 

226 

1260 

880.0 

20.99 

40.0 

1.6 

0.00 

226 

1290 

880.0 

19.71 

47.0 

1.5 

0.00 

226 

1320 

880.0 

20.07 

45.0 

1.5 

0.00 

226 

1350 

880.0 

19.29 

43.0 

1.4 

0.00 

226 

1380 

880.0 

18.12 

49.0 

0.7 

0.00 

226 

1410 

880.0 

16.39 

62.0 

0.7 

0.00 

226 

1440 

880.0 

17.04 

64.0 

1.0 

0.00 

Figure  4.  Sample 

data 

set. 
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TABLE  8.  COMPARISON  OF  MODELED  SENSIBLE  HEATS 
WITH  LITERATURE  DERIVED  VALUES 
(iinita  of  W/m*) 


Midday 

Night 

Surface  Type 

Model 

Literature 

Model 

Literature 

Snow 

75-125 

100 

-15 

-5 

Snow,  warm  advection 

65-100 

-50 

-13 

-100 

Wisconsin  farm 

170-260 

130 

-31 

-30 

California  farm 

169-270 

165 

not  avail. 

not  avail. 

Minnesota  farm 

120-210 

180-220 

not  avail. 

not  avail. 

Forest 

170-270 

175 

-34 

-49 

Desert 

250-350 

252 

-32 

-35 

Clay  pasture 

180-280 

215 

not  avail. 

not  avail. 

twenty  years.  The  results  used  in  the  modd  are  those  agreed  upon  by  most  investigators, 
especially  the  daytime  form  fit  that  was  thorouj^y  investigated  in  the  Patilson  (1970) 
paper.  The  nighttime  coefficients  are  much  less  universally  chosen,  and  though  a  P  value 
of  5  was  used  in  the  model,  a  change  to  7  or  8  may  be  warranted  (see  appendix  B). 

In  Tofsted  and  Gillespie  (1986),  TCHAD  (original  model)  was  compared  with  TCHAD  2.0, 
with  the  original  Deardorff  model  (here  named  SUHFA),  and  with  measured  data.  Figures 
5  through  10  are  reproduced  from  that  peq>er.  In  Figure  5  the  improvement  of  both 
TCHAD  and  TCHAD  2.0  over  SUHFA  is  seen  through  a  comparison  of  model  results 
of  sensible  heat  flux.  This  figure  shows  that  instituting  the  change  in  near-surface  air 
temperature  calculation  suggested  by  Kunkel  results  in  improved  sensible  heat  calculations. 
Both  the  original  TCHAD  and  TCHAD  2.0  are  superior  to  SUHFA.  SUHFA  results  in  an 
overestimate  of  sensible  heat,  while  simultaneously  restricting  heat  flow  in  other  forms. 
The  siirface  temperature  therefore  remains  high  into  the  nighttime  and  results  in  a  positive 
heat  flow  at  night  when  the  stuface  temperature  should  be  lower  than  the  air  temperature. 
Figure  6  shows  that  the  net  result  of  the  incorrect  heat  fliix  calculations  of  SUHFA  result 
in  large  errors  in  the  temperature  gradient  estimation.  Again,  TCHAD  and  TCHAD  2.0 
work  nearly  identically,  with  the  improved  version  performing  slightly  better. 

In  Figures  7  and  8  the  model  results  are  compared  with  a  time  history  of  the  refractive 
index  structure  parameter.  Figme  7  shows  the  thermal  probe  data,  which  normally  reads 
differently  than  scintillometer  data.  Figure  8  is  based  on  scintillometer  data,  and  modeled 
data  match  it  better.  Interestingly,  though  SUHFA  does  poorly  in  modeling  sensible  heat 
flux  and  temperature  gradient,  for  the  Deming  data  it  estimates  the  turbulence  data  as 
well  as,  if  not  better  than,  the  other  two  models  estimate  the  data.  But  when  compared 
to  the  White  Sands  data,  the  SUHFA  results  are  obviously  high  during  the  daytime  and 
in  the  early  evening  period. 

Figures  9  and  10  show  results  over  a  vegetated  field  (in  figure  9)  and  the  same  field 
when  snow-covered  in  winter  (figiure  10).  Again  SUHFA 's  ability  to  estimate  temperature 
gradient  (figure  9)  is  below  that  of  either  version  of  TCHAD.  In  figure  10,  SUHFA 


Figure  5.  Modeled  sensible  heats  are  compared  against  data  measured  at  Deming,  NM, 
on  14  Aug  1985  over  a  15  percent  foliage  cover  of  a  vineyard. 


Figure  6.  Measured  and  modeled  vertical  temperatiure  gradients  for  the  Deming  site 
14  Aug  85. 
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Figure  9.  Comparison  of  measured  and  modeled  temperature  gradients  for  a  site  near 
Flatteville,  IL,  on  24  Jun  1984. 
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Figure  10.  Comparison  of  measured  and  modeled  temperature  gradients  for  the  Flat¬ 
teville  site  on  17  Feb  1985. 
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overestimates  the  gradient  during  the  evening  period  because  low  windspeeds  unbalance 
the  equation  set  used  to  calcidate  the  flux-profile  parameters.  However,  TGRAD  does 
a  good  job  in  approximating  the  svimmertime  case,  but  fails  in  the  winter  night  case  — 
yet  for  different  reasons.  First,  TGRAD  does  not  use  the  near-surface  air  temperature 
method  of  TGRAD  2.0.  Second,  the  version  of  TGRAD  2.0  used  in  computing  this  figure 
allowed  the  characteristics  of  the  soil  to  vary  as  the  water  froze.  With  a  lower  thermal 
conductivity,  the  temperature  could  drop  faster,  and  the  gradients  increased  in  magnitude. 
But  even  TGRAD  2.0  was  incapable  of  predicting  the  correct  gradients  at  night,  due  to 
the  effects  of  condition  (67)  and  the  corresponding  need  to  set  u*  to  a  minimiim  value  to 
keep  the  algorithm  from  yielding  negative  L  values  at  night. 

The  case  of  high  temperature  gradients  at  night  is  equivalent  to  reaching  the  critical 
Richardson  number  where  turbulent  motions  in  the  atmosphere  are  no  longer  possible. 
There  are  several  possibilities  to  dealing  with  the  super-critical-Richardson-number  stable 
case.  One  method  tried  has  been  to  scale  the  gradient  upwards  based  on  a  functional  form 
related  to  the  amount  by  which  the  stability  exceeds  the  critical  value.  Another  method 
explored  would  be  to  find  a  stochastic  set  of  values  of  the  temperature  gradient  at  night. 
Therefore,  though  not  predicting  an  answer,  at  least  a  better  estimate  of  actual  gradients 
can  be  obtained  for  system  performance  purposes.  A  third  approach  would  be  similar  to 
the  second,  but  would  attempt  to  catagorize  the  statistics  by  timed  wave  actions  within 
the  surface  layer.  A  simulation  would  therefore  entail  the  setting  up  of  variotis  wave  states 
that  have  characteristic  break  times.  Each  wave  would  carry  a  certain  amoimt  of  the  total 
momentum  energy  of  the  surface  layer,  which  would  be  transfered  to  the  ground  in  the 
form  of  turbulence  at  the  time  of  the  wave  break.  The  net  result  would  be  a  combination 
of  empirical  observations  of  wave-structure-related-temperature-gradient-variations  and  a 
wave  model  that  would  start  after  sunset.  Unfortimately,  these  techniques  would  not  be 
sufficient  to  characterize  the  other  meteorolo^cal  parameters  that  must  be  modeled  to 
understand  the  surface  layer  atmosphere,  but  once  methods  had  been  developed  for  one 
variable,  similar  approaches  could  be  employed  for  the  others. 

6.  CONCLUSIONS 

The  model  TGRAD  has  been  shown  to  provide  accurate  estimates  of  daytime  surface 
layer  structures  for  temperature  gradient,  and,  by  implication,  for  other  surface  layer 
parameters  as  well.  Baseline  estimates  are  available  as  functions  of  the  time  of  day, 
including  nocturnal  temperature  structure  estimates.  The  model  derived  is  a  significant 
improvement  in  performance  over  the  previous  Deardorff  model  on  which  it  is  based,  while 
retaining  the  feature  of  estimating  sensible  heat  flux,  taking  into  account  foliage  effects. 
The  model  is  capable  of  estimating  all  the  flux-profile  parameters  necessary  to  calculate  the 
refractive  index  structure  parameter  as  well  as  the  vertical  temperature  and  wind  profiles. 

At  night,  the  problem  of  a  highly  stable  atmosphere  is  significant,  and  as  yet  not  completely 
solved.  Methods  have  been  described  whereby  progress  in  this  area  may  be  made.  For  now 
the  model  described  is  sufficient  to  correctly  estimate  temperature  gradients  up  to  about 
0.5  *’C/m  at  2  m  a.s.l. 
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APPENDIX  A 

CALCULATIONS  INVOLVING  SATURATION  SPECIFIC  HUMIDITY 


The  specific  htunidity  at  a  given  temperature  b  fimnd  by  multiplying  the  saturatkm  specific 
humidity  {Q»mt)  hy  a  fraction  between  zero  and  one,  where  thb  fraction  will  either  be  the 
relative  humidity  or  the  moisttue  coefficients  of  the  ground  or  vegetation. 

Q,mt  b  normally  written  as  a  function  of  pressure,  P,  and  saturation  vapor  pressure,  e««t. 
The  saturation  vapor  pressure  is  a  function  of  the  temperature. 

=  6.1078  X  (A  -  1) 

where  over  water  a  is  7.5  and  6  is  237.3,  while  over  ice  a  is  9.5  and  b  is  265.5.*  Q,Mt  b 
rebted  to  Cgat  and  P  through  the  equation 


Q»at  — 


0.62  Cmi 

P- 0.38  c..,’ 


(A -2) 


which  is  also  found  in  Haurwitz.* 

The  derivative  of  Q..,  with  respect  to  temperature  is  a  function  of  the  derivative  of  e.., 
with  respect  to  temperature: 


Thus 

'W  "  (6+r)*  { i-o.38c,.,/p  }  ~ 

This  result  was  found  through  taking  the  derivative  of  Q..,  in  terms  of  de..(/dT  and 
finding  de.a,/dT  through  taking  the  derivative  of  Inc..,  =  ln(10)/(T)  and  solving  for 
dc.«,/<ir . 

Pressure  is  treated  as  constant  with  temperature  since  it  b  an  input  dependent  on  synoptic 
events  rather  than  on  the  temperatures  treated  within  the  model.  The  derivative  of  Q.., 
with  respect  to  temperature  is  thus  only  a  function  of  e..,. 


*Haurwitz,  B.,  1941,  Dynamic  Meteorology,  McGraw  Hill,  New  York,  pp.  8-10. 
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APPENDIX  B 


CONSIDERATION  OF 

NOCTURNAL  TEMPERATURE  GRADIENT  STRUCTURES 


1.  OVERVIEW 

Surface  energy  budget  methodologies,  combined  with  flux-profile  equations,  can  be  used 
to  estimate  the  vertical  temperature  structure  of  the  surface  boundary  layer  from  data 
observed  at  a  single  height  during  the  day.  However,  at  night,  the  governing  equations 
in  this  method  often  fail  as  the  critical  Richardson  number  is  surpassed  (Hansen,  1977). 
Under  these  conditions,  gravity  waves  propagating  within  the  inversion  layer  are  crucial  in 
understanding  temporal  fluctuations  of  temperature  gradient  strength. 

In  this  appendix,  a  model  is  proposed  to  characterize  the  break  frequencies  of  these  wave 
modes.  In  this  model  waves  are  considered  to  reflect  at  both  the  surface  and  close  to  the  top 
of  the  inversion.  Also,  the  vertical  gradient  of  horizontal  wrind  and  the  second  derivative 
of  temperature  are  treated  as  being  constant  with  height.  The  vertical  wavelength  is 
considered  to  be  an  integral  fraction  of  the  inversion  height,  and  the  horizontal  wavelength 
is  considered  to  be  7.5  times  the  depth  of  the  inversion  layer. 

Model  predictions  of  frequency  modes  at  0.210,  0.109,  and  0.074  cycles  per  xninute 
(cyc/min),  and  corresponding  break  frequencies  of  0.128,  0.034,  and  0.011  cyc/min  appear 
to  be  closely  matched  to  test  data  values. 

2.  INTRODUCTION 

Over  the  years  many  studies  have  treated  the  problem  of  strong  stability  conditions  within 
the  surface  layer  atmosphere.  But  the  noct\imal  problem  has  continued  to  present 
theoretical  difiiculties.  In  particular,  the  application  of  flux  profile  parameterizations 
does  not  hold  during  highly  stable  periods.  These  parameterizations  depend  on  a  fully 
turbulent  atmosphere  for  their  underlying  theoretical  basis.  But  at  night,  calm  periods 
occur  when  turbulent  mixing  is  highly  damped. 

In  the  absence  of  the  flux-profile  parameterizations  of  stable  vertical  atmospheric  structure, 
some  other  model  must  be  postulated.  As  evidenced  by  equation  (61)  in  the  main  section  of 
this  report,  the  flux-profile  technique  caimot  obtain  a  valid  solution  when  the  temperature 
inversion  strength  impedes  vertical  mixing  of  air. 

As  was  discussed,  the  quick  fix  to  the  problem  was  to  limit  the  computation  of  u*  to  some 
minimum  value.  In  this  way  the  resulting  algorithm  always  obtains  a  positive  value  for  the 
Obukhov  length.  However,  this  does  not  ensure  a  valid  estimation  of  the  vertical  structure 
of  the  temperature  profile.  Perhaps  this  is  because  the  vertical  temperature  profile  is 
relatively  independent  of  the  flux-profile  forms  under  strong  inversion  conditions.  This  is 
the  thesis  of  this  section. 


When  the  predominant  processes  governing  the  vertical  structures  of  the  main  parameters 
within  the  surface  layer  are  not  due  to  wind  driven  mixing,  there  is  no  readily  available 
theory  to  gain  predictive  insight  into  these  structures.  An  alternative  driving  mechanism 
must  therefore  be  sought.  It  is  hypothesized  that  there  are  two  driving  mechanisms:  the 
radiative  flux  at  the  surface  that  is  cooling  the  near  surface  atmosphere,  and  standing 
gravity  waves  within  the  boundary  layer.  The  gravity  waves  so  produced  tend  to  gain 
energy  exponentially  due  to  the  nonlinear  nature  of  the  governing  wave  equations.  Then, 
since  each  wave  is  subject  to  a  turbulent  breakdown  and  since  this  breakdown  has  a 
characteristic  time  for  each  wave,  the  dynamic  variability  of  vertical  temperatvue  gradients 
with  time  may  be  explained  by  periodic  mixing  episodes  corresponding  to  gravity  wave 
breakdown. 

It  is  hoped  that  once  the  temporal  structure  of  the  vertical  temperatme  gradient  is 
understood,  then  the  results  can  be  coupled  to  models  indicating  the  time  variability 
of  the  flow  of  turbulent  energy  to  the  surface  and  the  dissipation.  The  knowledge  of  the 
energy  flow  should  help  in  the  calculation  of  the  parameters  u,,  0*,  and  L. 


3.  ENERGY  BALANCE  CONSIDERATIONS 

Following  the  traditional  approach  to  the  computation  of  u*  and  L,  the  nocttunal  energy 
budget  is  divided  into  components  related  to  longwave  radiative,  conductive,  and  convec¬ 
tive  fluxes.  Since  in  this  analysis  we  are  considering  only  the  most  stable  atmospheric 
states,  we  can  safely  assume  the  radiative  and  conductive  fluxes  will  be  the  dominant 
drivers  in  the  energy  budget.  Rewriting  equation  (3), 

+  Hsg  +  LEg  +  G  =  Rifi  {R  ~ 

In  this  model  the  only  means  of  cooling  the  air  near  the  surface  is  through  the  sensible  heat 
flux,  but  under  quiescent  conditions  it  may  be  that  radiative  exchange  of  energy  with  the 
surface  will  cause  cooling  of  the  air  near  the  surface.  This  cooling  will  depend  on  numerous 
factors  beyond  the  scope  of  this  section.  However,  these  additional  factors  may  play  key 
roles  in  the  stable  boundary  layer  problem.  Such  factors  include,  but  may  not  be  limited 
to,  radiative  flux  divergence  in  the  stable  layer,  mechanisms  for  generating  the  nocturnal 
jet  flow,  drainage  winds,  and  vertical  stratification  of  moisture  in  the  stable  layer. 


4.  THE  CRITICAL  RICHARDSON  NUMBER 
At  night  Dyer  (1974)  finds 

V'fc  =  V’m  =  V’ = ,  (^-2) 

as  pointed  out  in  the  main  text.  However,  Hansen  (1977)  also  indicates  the  choice  of  the 
critical  value  of  the  gradient  Richardson  number  is  related  to  ^  through 

0  =  I//?....,  .  (B  -  3) 
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This  critical  value  is  related  to  the  value  beyond  which  turbulent  mixing  is  thought  to  be 
damped.  Note  that  this  condition  is  the  finding  of  the  analysis  in  the  main  text  leading 
to  condition  (67).  Accordingly, 

Ri  =  gzie,  - 


firom  the  analysis  in  the  main  text,  since  when  Ri  reaches  its  critical  value  condition  (67) 
will  be  met.  Another  way  of  expressing  this  equation  is 


Ri  = 


gT»  z 

U«  U$x' 


where  T*  is  the  scaling  temperature  (Hoffert  and  Storch,  1979),  given  by 


^  k{9x-Bx) 

*  MC)  -  0  ’ 

and  u«  is  the  friction  velocity  given  by 

_ 

€n(C)  - ' 

Assuming  T«  and  u«  can  both  be  determined,  according  to  the  flux-profile  theory  they 
are  considered  constant  over  the  surface  layer.  Due  to  this  relationship,  Ri  is  a  function 
of  height.  Also  we  can  see  that  since  u  is  zero  hi  z  =  zq,  Ri  is  infinite  at  zero  height 
and  decreases  thereafter.  We  could  therefore  expect,  even  under  the  limitations  of  the 
flux-profile  method,  that  there  is  some  height  above  which  mechanical  turbulence  exists. 
The  question  that  remains  is  how  to  treat  the  intervening  layer. 

In  answering  this  question  the  critical  Richardson  number  is  a  valuable  meastire  of  the 
stability  of  the  stable  nocturnal  boundary  layer,  and  a  study  of  turbulence  formation  in 
this  layer  naturally  arrives  at  an  estimation  of  this  parameter.  Such  a  study  is  performed 
in  the  next  section  with  the  help  of  a  mathematical  construct  called  a  turbulence  rotor. 


5.  A  TURBULENCE  ROTOR 

We  can  gain  further  insight  into  the  cutoff  of  convective  fluxes  and  determine  a  theoretical 
value  for  the  critical  gradient  Richardson  number  if  we  model  a  single  rotor  of  turbulence 
within  the  flow.  Consider  a  portion  of  the  general  flow  of  air  above  the  surface.  Within 
a  certain  small  height  interval,  the  flow  can  be  described  by  average  vertical  gradients  of 
windspeed  and  temperature.  It  is  postulated  that  a  parcel  of  air  cjw  become  turbulent 
if  it  becomes  disconnected  from  the  general  flow  and  begins  to  roll  over  (assiuning  this 
disconnection  is  through  the  addition  of  minimal  startup  energy).  Assuming  random 
small  scale  turbulent  energy  exists  throughout  the  flow,  such  a  parcel  is  always  able  to 
arise  spontaneously  from  within  the  flow. 
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A  further  postulation  is  that  the  simplest  such  starter  parcel  would  be  shaped  like  a 
cylinder.*  Its  time  of  existence  would  be  short:  drag  forces  drawing  it  back  into  the  flow 
in  time.  But  during  its  existence  it  would  decrease  the  vertical  gradients  of  temperature 
and  momentum  through  mixing  of  the  air  at  different  temperattures  and  velocities  at  its 
upper  and  lower  halves.  Figure  B-1  shows  the  model  for  the  hypothesized  cylinder  at  the 
time  of  its  detachment  from  the  flow. 


Figure  B-1. 


2  axis  (height) 


Vertical  cross  sections  of  temperature  and  wind  relative  to  a  turbulence 
cylinder. 


At  disconnect  time,  two  competing  factors  are  influencing  the  cylinder:  the  rotational 
kinetic  energy  possessed  by  the  cylinder  that  causes  it  to  roll  and  the  negative  potent’ 
energy  that  tends  to  prevent  it  from  rolling.  The  first  of  these  may  be  defined  by  computing 
L,  the  cylinder’s  angular  momentum,  as 


«2, 


(B-4) 


where  M  is  the  mass  of  the  cylinder,  dm  is  a  differential  mass  element,  r|,  =  r  8in(n7)  is  the 
vertical  distance  away  from  the  center  of  the  cylinder,  and  u,  is  the  horizontal  windspeed 
at  height  r|x.  (See  figure  B-2  for  r  and  w  definitions.) 

Then,  with  po  the  air  density  at  the  height  of  the  center  of  the  cylinder,  uq  the  windspeed 
at  that  height,  and  assuming  p  and  u  (the  vertical  derivatives  of  density  and  horizontal 
windspeed  across  the  cylinder)  are  constant,  we  have 


«2  =  tto  +  «r8in(tn), 


*Such  structures  were  observed  in  experiments  where  a  collimated  HeNe  laser  was  propa¬ 
gated  through  a  2000  m  path  of  air  on  a  calm  night.  In  most  cases  these  turbulent  cylinders 
or  spheres  appeared  to  ‘roll’  through  the  illuminated  surface  where  the  laser  pattern  was 
intercepted. 
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Figure  B*2.  Definitions  of  radius  (r),  an^e  (tn),  and  width  (y)  of  a  turbulence  cylinder. 


dm  =  pdV  prdrdadyy 
P  =  P(t  +  prsm{w), 
p  =  dpfdz, 
li  =  du/dz. 

These  approximations  become  exact  as  the  vertical  size  of  the  tylinder  approaches  zero. 
Also,  there  should  be  no  confusion  between  this  notation  and  the  standard  use  of  the  over¬ 
dot  as  a  time  derivative  since  time  derivatives  are  not  explicitly  made  in  this  derivaticm. 

Equation  (B-4)  can  now  be  rewritten 


V  2ir  R 

=  ///{(^ 


+  p  r  sin(m))  r  dr  dw  dy){r  8in(to)}  {«o  +  «  r  sm(tn)} 


0  0  0 


=  (po«  +  puo)—  ,  {B-b) 

where  V  is  the  volume  of  the  cylinder.  The  three  quantities  within  brackets  inside  the 
integral  correspond  to  the  same  three  quantities  in  equation  (B-4). 

From  the  angular  momentmn,  the  rotational  kinetic  energy  can  be  calculated. 


where  I  is  the  rotational  inertia.  For  a  cylinder, 

j 

2 


73 


so  that  2 

£k  =  |{W»  +  ^»o)^}  /(>.Vr*  .  {B-6) 

The  next  quantity  to  be  determined  is  the  rotational  potential  energy  to  be  overcome.  To 
compute  this  quantity  the  torque  must  be  evaluated.  The  tOTque  induced  by  rotating  the 
cylinder  in  the  presence  of  a  density  gradient  can  be  defined  relative  to  an  arbitrary  angle 
the  angle  of  rotation  from  the  equilibrium  position.  Computing  the  torque  involves 
another  integration  over  the  cylinder’s  voltune. 

(B-7) 

Based  on  the  geometry  of  figure  B-2,  the  torque  can  be  defined  as  t  J  =  r,  where  j  is  the 
unit  vector  in  the  y  direction.  In  simplifying  the  notation,  the  torque’s  y-component  will 
hereafter  be  referred  to  as  the  torque,  r.  The  right-hand  side  of  equation  (B-7)  can  be 
calculated  using  the  figure  B-2  definition  of  r,  and 


f=/, 


xdF. 


dF  =  —dpgdV  ib, 


where  dp  is  the  density  difference  between  the  volume  element’s  equilibrium  level  and  its 
current  level,  and  k  is  the  unit  vector  in  the  +z  direction.  Since  density  is  assumed  to 
v’ary  linearly  with  height, 


dp  =  p  r  (sin(m  -ft?)  —  sin(t!7)). 


y  2ir  R 


-III  r  cos(c7  f  1?)  {—pr  (sin(m  ft?)  —  sin(m))  grdr  dw  dy) 


0  0  0 


=  V’sin(t?) 


(B-8) 


Assume  that  the  cylinder  can  be  considered  a  true  turbule  if  it  can  rotate  through  tr  radians 
(it  can  turn  over).  The  amount  of  potential  energy  that  must  be  overcome  is  fotmd  by 
integrating  the  torque  through  tt  radians. 


Ep  = 


/rdt?  =  2pyV-^  . 

0 


Setting  the  amoimt  of  potential  energy  to  be  overcome  equal  to  the  kinetic  energy  available, 
an  equation  for  the  necesszoy  tit  is  developed. 


0  =  u*  -I- 


8£p\ 
fio  ) 


(B-9) 
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Solving  for  u, 


(B  - 10) 


u  =  ^6gQ~  uoQ, 

where 

0=f.  (B-ll) 

po 

FVom  this  development,  u  defines  the  minimum  positive  gradient  of  windspeed  necessary 
for  turbulence  to  be  initiated  (ignoring  drag  effects  that  might  be  present  over  the  first  sr 
radians).  Farther,  since 

Po  Bdz* 

(ignoring  vertical  pressure  gradient  effects  in  a  highly  stable  atmosphere),  therefore. 


du  _  ISg^  uM 
dz  y  B  dz^  Bdz 


{B  - 12) 


where  g  =  9.8  m/s*. 

To  simplify  the  above  expression,  consider  the  following  example:  In  a  strongly  stable 
atmosphere  let  z  =  2  m,  ^  =  290  K,  u  =  2  m/s,  and  dBfdz  =  0.5  K/m.  Then, 

~  =  0.368s“*  +  0.00345s’*. 
dz 

Clearly,  the  first  term  is  the  most  important,  and  if  the  second  term  also  becomes  significant 
because  of  higher  windspeeds,  the  flux-profile  methods  will  not  fail.  Therefore  only  the 
first  term  needs  to  be  retained. 

Therefore,  keeping  only  the  first  term  on  the  right,  squaring  what  remains  of  equation 
(B-12),  and  moving  terms  result  in 


9^  / ( duV 
Bdz/\dz)\  ’ 

t 


(B  - 13) 


where  the  >  sign  is  introduced  to  indicate  that  the  condition  b  necessary  for  turbulence. 
But  the  term  on  the  right  is  the  same  equation  given  by  Hansen  for  the  gradient  Richardson 
number,  implying  a  critical  value  under  very  stable  conditions  of  j,  and  thus  $  equals  8. 

This  result  corresponds  well  with  the  conclusions  of  Oke  (1970)  and  McVehil  (1964) 
that  turbulence  drops  sharply  aroimd  Richardson  numbers  of  0.1  and  0.08,  respectively. 
Since  0.125  is  an  upper  limit  according  to  the  theory  developed  here,  test  data  at  lower 
stabilities  will  indicate  lower  effective  critical  values  (through  equation  (B-10))  in  accord 
with  Oke  and  McVehil.  But  this  result  does  not  compare  well  with  the  conclunons  of 
other  researchers  such  as  Webb  (1970)  and  Businger  et  td.  (1971),  who  have  found  critical 
Richardson  numbers  around  0.20,  though  there  is  no  complete  agreement  cm  this  subject 
(see  particularly  Yaglom  (1977)  and  Hansen  (1977)). 
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To  imderstand  the  discrepancy  between  equation  (B-13)  and  these  other  findings,  first 
consider  that  the  arguments  addressed  thus  far  have  d^Jt  only  with  turbulence  arising 
from  the  mean  fiow  characteristics.  Since  turbulence  is  conuncmly  observed  at  higher 
Richardson  numbers,  other  sources  of  turbulence  must  exist  that  provide  turbulent  energy 
between  the  general  fiow  cutoff  value  commonly  understood  to  be  near  0.2. 

One  such  source  of  turbulence  appears  to  come  from  wind  interaction  with  roughness 
elements  at  the  surface.  We  can,  in  fact,  argue  for  such  surface  based  turbulence  on 
similarity  theory  grounds. 

Assuming  that  the  windspeed  and  temperature  profiles  are  similar  (as  similarity  theory 
maintains),  the  ratio  of  the  gradients  of  temperature  and  wind  should  be  equal  regardless  of 
height  above  the  surface  at  night.  According  to  Dyer,  similarity  theory  expresses  nocturnal 
gradients  of  wind  and  temperature  as 


(B  -  14) 

(B  - 15) 

Therefore,  equation  (B-13)  can  be  reexpressed  as 

gkT.  LJ^ 

(B  - 16) 

Assuming  that  meaningful  values  of  d*,  u„  and  L  can  be  drawn  from  profile  data,  there 
must  be  some  level,  z,  low  enough  that  the  Ijz  term  in  equation  (B-16)  is  dominant,  and 
the  inequality  can  be  satisfied.  Turbulence  elements  produced  between  the  zero  plane  and 
this  level  will  then  generate  turbulence  that  can  propagate  upward  into  the  turbulence 
suppressed  region.  This  effect,  however,  will  only  influence  near-surface  conditions,  which 
would  not  allow  for  transfer  of  energy  between  the  base  and  top  of  the  inversion  layer. 

A  second  source  of  continuing  turbulence  may  be  remaining  energy  from  the  unstable 
boundary  layer  in  existence  before  sunset.  As  noted  by  Arya  (1972),  preexisting  turbulence 
must  be  considered  when  turbulence  at  Richardson  numbers  of  0.25  or  more  is  observed. 
However,  even  this  turbulence  will  rise  to  the  top  of  the  stable  layer  with  time,  and  thus 
does  not  provide  a  mechanism  for  maintaining  turbulent  energy  exchange  with  the  siirface 
throughout  the  night. 

6.  GRAVITY  WAVES 

A  third  source  for  turbulence  is  caused  by  the  breakdown  of  waves,  which  allows  transport 
of  tiurbulent  energy  to  the  surface  throughout  the  stable  period.  These  waves  are 
characterized  by  periodic  breakdowns  that  distmrb  the  surface  layer  atmosphere. 

Kelvin- Helmholtz  waves  (Atlas,  1970)  and  other  types  of  gravity  waves  (Davis,  1976) 
produced  throughout  the  inversion  layer  grow  and  periodically  break,  producing  turbulent 
episodes  diiring  an  otherwise  calm  evening.  These  periods  often  occur  as  dips  or  spikes  of 
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Figure  B-3.  Vertical  temperature  gradient  at  2  m  above  grotmd  level  as  a  function  of  time. 

This  is  period  2  of  data  taken  at  WSMR,  NM,  on  the  night  of  8  Sep  1984. 

ttirbulence  throughout  an  evening  (for  example,  see  Kunkel  and  Walters’  (1983)  figures  1 
through  3).  Similarly,  the  time  structure  of  the  temperature  gradient  does  not  vary 
smoothly  during  the  night,  as  might  otherwise  be  implied  by  Wyngaard’s  study  (1975). 
These  waves  are  manifested  in  periodic  increases  and  decreases  in  the  vertical  temperature 
gradient  (figure  B-3). 

Some  of  the  key  findings  from  studies  conducted  to  understand  gravity  wave  behavior 
are  presented  in  this  section.  A  general  treatise  on  gravity  waves  is  available  in  Beer 
(1974),  where  the  elemental  energy  and  momentum  conservation  equations  are  used  to 
derive  several  wave  equations  under  different  atmospheric  conditions.  A  modified  version 
of  Beer’s  method  for  an  inversion  layer  case  will  be  developed  here.  Beer  notes  that  gravity 
waves  are  reflected  at  a  level  where  their  fundamental  firequency  matches  the  Brunt- Viasala 
frequency.  This  means  that  all  waves  propagating  within  the  surface  layer  will  remain 
within  this  layer,  reflecting  at  the  ground  as  as  near  the  top  of  the  inversion. 

However,  a  fraction  of  these  waves  is  never  reflected  since  they  are  absorbed  before  reaching 
the  inversion  height.  The  height  of  absorption  is  characterized  by  the  condition  (Booker, 
1967) 

ui  —  KgU»  =  0  ,  (B  — 17) 

where  u  is  the  radial  wave  frequency.  Kg  is  the  horizontal  wave  ntunber,  and  u*  is  the 
windspeed  at  the  absorption  height.  The  left-hand  side  of  the  equation  is  also  called 
the  Doppler  velocity,  D.  Thus,  absorption  occurs  at  the  level  where  the  Doppler  velocity 
vanishes.  This  covers  most  short  wavelength  modes  (large  wave  number)  that  are  absorbed 
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well  below  the  inversion  level.  For  longer  wavelength  modes,  there  is  the  possibility  for 
tunneling  through  the  potential  barrier  at  the  inversion  height,  but  this  case  b  not  dealt 
with.  Results  obtained  by  ignoring  tunneling  are  virtually  identical  to  a  full  analysb,  but 
are  mathematically  much  simpler  (Orlanski,  1973). 

A  rich  selection  of  material  is  avmlable  on  the  rate  of  growth  of  different  wavelength  modes. 
Results  vary,  depending  on  the  density  and  windspeed  vertical  structures  chosen,  but  as 
Lcdas  and  Einaudi  (1976)  report,  the  wavelength  of  the  fastest  growing  mode  has  been 
fotmd  as  being  between  7.14  and  7.66  times  the  depth  of  the  shear  layer  studied.  Davis 
and  Peltier  (1976)  find  a  fastest  growing  wavelength  of  6.55  times  the  depth  of  the  inversion 
layer.  Traditionally,  however,  the  coefEcient  7.5  appears  most  frequently  in  the  literature 
(Browning,  1971;  Miles  and  Howard,  1964;  Goldstein,  1931),  and  apparently  dates  back  to 
the  works  of  Rayleigh  (1880, 1887).  This  traditional  value  was  used  in  this  analysis.  In  our 
case  the  characteristic  shear  layer  depth  is  the  depth  of  the  inversion  layer,  Zi,  assumed 
to  range  between  100  and  400  m  (Arya,  1981).  Corresponding  fastest  growing  mode  wave 
numbers  are 

■‘*  =  7^-  • 


6.1  Derivation  of  the  Wave  Equation 

The  derivation  of  the  wave  equation  governing  surface  layer  propagating  gravity  waves 
closely  follows  a  procedure  described  by  Beer  (1974).  In  his  exposition.  Beer  manipulated 
the  equations  of  motion  in  x  and  z  directions,  x  being  along  the  direction  of  the  wind 
and  z  being  the  upward  direction.  He  also  included  the  energy  conservation  and  mass 
continuity  equations,  to  obtain  four  equations  similar  to  those  of  equation  set  (B-19). 


-iK^gH  (g)  + 

“'.^+  Q<  =  o 

(B-19a) 

+  5u'.  =  0 

(B-19b) 

7X- 

^  +  §(g)=o 

(B-19c) 

X  - 

4(i  +  ^)  +  §(^)=0  , 

(B-19d) 

where  i  =  \/^  is  the  imaginary  root,  H  is  the  atmospheric  density  scale  height  (H  w  8000 
m),  X  is  the  wind  divergence, 

du'  du' 

^  =  ■ 

Q  is  &  linearized  operator  pro])ortional  to  the  Doppler  frequency  {Q  =  tfi),  7  is  the  ratio 
of  specific  heats,  and  the  four  variables:  pressure  (P),  density  (p),  horizontal  windspeed 
(u*),  and  vertical  windspeed  (u,),  are  considered  to  consist  of  an  average  value  and  a  small 
perturbation  (the  primed  portion).  For  example 


Equation  set  (B-19)  differs  from  Beer’s  set  (equations  (2.5.1)  through  (2.5.4))  by  the  last 
term  on  the  left-hand  side  of  equation  (B-19a).  This  tenn  accounts  for  the  gradient 
of  horizontal  windspeed  with  height  in  the  surface  layer.  Beer’s  analysis  results  in  the 
following  wave  equation 


dz^ 


-M' 


(K2c2-Q2) 


dH\  dW 
\  dz  }  dz 


K*c* 


7fr*(Kic2-ft2) 


By  using  the  following  substitutions 

'fgH  =  c'^ 

0  =  KI<? 

equation  (B-21)  can  be  rewritten  as 

0 


(B  -  21) 


ii  ^  iH 

"  =  7? 


+  +  =  0  -  (B-22) 


6.2  Wave  Frequencies  and  Break  Intervals 

W  represents  that  portion  of  the  solution  for  the  vertical  perturbation  velocity  that  is 
dependent  on  height.  In  Beer’s  approach  the  proposed  form  for  u',  is 

u',(x,z,t)  =  W(z)  exp(i(u;<  -  Kxx))  ,  (B-23a) 

where 

W{z)  =  Z  exp(-t  K.  z)  (B  -  236) 

=  >  (B  -  23c) 

and  Z  is  a  complex  constant.  W  therefore  represents  the  envelope  of  the  complex  phase. 
This  approach  assumes  the  equation  of  motion  is  separable.  Also,  it  assumes  fixed  forms 
for  the  vertical  temperature  and  mean  windspeed  structures. 

To  find  a  frequency  u;  for  this  form  of  solution,  equation  (B-22)  can  be  rewritten  as 

W  +  DW  +  SW  =  {i  .  (B-24) 

Then  using  equation  set  (B-23),  the  wave  equation  can  be  broken  into  two  equations, 
written  in  terms  of  real  and  imaginary  parts  of  the  wave  number  K^. 

+£>«,. +5  =  0  (B-25a) 
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{B  -  256) 


2*(k,,  =0 

Solving  these  equations  for  k,.  and  k,,,  we  find 

D 


pa  . 


By  using  the  boundary  conditions 


(B-26a) 

(B-266) 


W^(0)  =  W^(Z0  =  0  , 

we  obtain 

=  ^  =  ys+jc*,  n  =  l,2,3...  ,  (.8-27) 

which  is  a  similar  result  to  that  obtained  by  Davis  and  Peltier  (1976),  where  standing 
waves  appear  below  the  inversion  level. 

This  discussion  can  be  extended  to  equation  set  (B-19)  quite  readily  once  a  wave  equation 
comparable  to  equation  (B-22)  is  found.  Working  throu^  equation  set  (B-19)  with  Beer’s 
method  results  in 


j  0  -  /c,  ti,  UH't)  ( 


r  ;  'tgH 


|w  = 


(B-28) 


By  inserting  the  new  forms  for  D  and  5  into  equation  (B-27),  a  transcendental  equation 
in  uj  emerges  that  can  be  solved  iteratively  with  the  Newton-Raphson  method  as  follows: 
Set 


u  = 


Prom  equation  (B-27)  a  function  J(u;)  can  be  used 


J(u;)  =  -I- 5  —  1/ =  0  , 

whose  roots  must  be  found.  These  roots  are  determined  by  selecting  some  initial  value  for 
u):  u^Q.  The  Newton-Raphson  method  is  then  used  to  determine  the  root: 


Wm+l  =  Wm  -  J{ujfn) 


/(£) 


(Wm) 
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where 


dJ  ZDdD  ^  dS 
dui  2  du  du 

—  -  “«> 

du  ^njy  ft» 


ds  ^  a 

du)  'YgH 


ra  jj  **  “**  7  “*v  ns  n'r*\2 


Using  this  method  the  value  of  ui  can  be  computed  at  a  given  height.  In  this  wave  equation, 
however,  the  value  of  u)  will  vary  with  hdght.  This  characteristic  of  the  wave  leads  to  an 
interesting  paradox:  while  Q  was  assumed  constant  with  height,  u;  is  a  function  of  height 
whose  gradient  is  not  necessarily  small.  The  resolution  to  this  apparent  contradiction  is 
that  when  the  governing  equations  are  reviewed,  we  find  that  the  terms  involving  dSl/dz 
are  small  compared  to  other  factors  in  the  equations.  Thus  the  wave  frequencies  computed 
as  functions  of  height  will  not  be  significantly  in  error. 


The  reason  for  a;’s  variation  with  height  comes  from  the  necessity  that  be  conserved 
throughout  the  depth  of  the  inversion  layer.  If  this  were  not  true,  the  wave  would 
destructively  interfere  with  itself  at  every  reflection  and  little  wave  energy  could  be  built 
up.  The  wave  frequency  therefore  changes  with  hdght  so  that  k,,  can  be  “transported” 
within  the  layer,  as  Whitham  (1974)  points  out  in  detail. 

Following  Whitham’s  arguments  further,  the  wave  propagation  established  when  any  wave 
has  a  frequency  that  varies  with  position  is  inherently  unstable.  That  is,  the  wave  will 
eventually  catastrophically  interfere  with  itself,  creating  a  brealdng  action  similar  to  that 
experienced  by  a  wave  on  the  ocean  as  it  approaches  the  shore.  The  crest  of  the  wave  rises 
precipitously,  and  cannot  support  itself.  The  natiu«  of  the  characteristic  time  of  existence 
of  such  a  wave  is  a  function  of  the  initial  differences  in  phase  between  the  top  and  the 
bottom  of  a  wave.  Assuming  the  wave  is  in  phase  at  the  time  of  its  building,  it  proceeds 
toward  a  breakdown  at  a  characteristic  time  later  that  Whitham  describes  as 


=  - 


(B-29) 


Computationally,  this  equation  can  be  determined  by  finding  u;  for  two  close  values  of 
around  the  solution  values  for  these  two  terms  at  some  hdght  z.  Then,  dufldKg^  can  be 
evaluated  at  a  second  height  z  +  dz  and  the  derivative  computed. 

When  this  procedure  was  followed  for  temperature  gradient  data  collected  at  WSMR  in 
1984  and  compared  with  the  Foiuier  transform  of  this  data,  it  was  discovered  (Tofsted, 
1987)  that  wave  frequencies  and/or  breakdown  frequendes  were  present  in  the  data  at 
aroimd  0.0018,  0.004,  0.008,  0.012,  and  0.015  radians/sec  (0.017,  0.038,  0.076,  0.115,  and 
0.143  cydes  per  minute).  As  mentioned  in  the  introduction,  modd  predictions  were  made 
to  detomiine  characteristic  wave  frequendes  for  this  case.  FVequendes  of  0.210,  0.109, 
and  0.074  cydes  per  minute  and  corresponding  break  frequendes  of  0.128,  0.034,  and 
0.011  cyc/min  were  determined  based  on  the  mean  temperature,  temperature  gradient, 
and  windspeed  observed  during  these  trials. 
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This  case  shows  that  apparently  there  is  a  reasonable  connection  between  the  observed 
data  and  the  predictions.  More  trials  would  be  necessary  to  validate  this  hypothesis.  Also, 
the  method  does  not  help  in  providing  any  predictive  infonnation  concerning  the  expected 
strength  of  the  gradient  and  the  interference  effects  between  energy  released  when  one 
wave  breaks  the  buildup  of  energy  at  another  wavelength. 

The  source  of  gravity  wave  energy  apparently  is  the  wind  field  aloft.  Lilly  and  Klemp  (1979) 
and  Peltier  and  Clark  (1979)  focus  on  mountain  wave  formations,  but  orographic  effects  of 
nearly  any  degree  are  probably  sufficient  to  produce  waves  effects  of  all  wavelengths  due 
to  the  dispersive  nature  of  gravity  waves.  In  particular,  for  regions  near  mountains,  lee 
waves  can  transport  energy  downward  into  the  inversion  layer  where  it  becomes  effectively 
trapped  through  conversion  into  shorter  wavelength  modes.  Once  confined,  these  waves 
are  capable  of  carrying  momentum  vertically  within  the  inversion  layer  without  turbvdent 
interaction.  As  noted  by  Arya  (1981),  this  wave  energy  could  be  the  source  of  continuing 
turbulence  up  to  Richardson  numbers  of  0.25,  since  at  the  time  of  wave  breaking,  energy 
carried  within  the  waves  should  convert  directly  into  localized  turbulence.  Also,  in  view  of 
dispersion,  there  should  always  be  background  turbulence  fin>m  energy  transfer  to  shorter 
wavelength  modes.  With  very  short  breaking  times,  these  modes  continually  feed  small 
amounts  of  turbulent  energy  into  the  flow. 

Using  the  gravity  wave  model  one  possibility  for  simulating  fluctuations  in  vertical  tem¬ 
perature  gradient  would  be  to  treat  the  wave  mode  effects  as  independent  and  allow  them 
to  influence  the  total  temperature  gradient  as  a  sum  of  effects.  This  technique  would  use 
an  equation  of  the  type 

^  sin(wit  +  e„,)F(lB„  08,,  t)  +  ^,  (B-  30) 

i 

where  t)  would  describe  function  for  the  bounding  envelope  of  the  building  and 

breaking  of  the  wave,  and  Ow,  and  0  b,.  are  random  phases  of  the  wave  and  breaking  times. 
Simple  functions  describing  the  F  function  would  be  sines  and  cosines  or  combinations 
thereof.  Similar  functions  could  be  developed  for  Cf,  C^,  and  C*.  Since  the  i>eriods  of 
breaking  of  some  of  these  waves  are  as  little  eis  5  min,  measurements  must  be  taken  at 
least  10  times  as  fast  as  this  rate  and  the  time  scale  of  the  surface  energy  budget  model 
must  be  set  as  fast  or  faster  than  this  rate.  Thus  at  night  the  time  step  should  be  no  more 
than  30  seconds  (when  wave  phenomena  are  considered). 


7.  Approximation  of  Temperature  Gradient 

In  lieu  of  a  temperature  gradient  calculation  based  on  gravity  wave  effects  (which  is 
currently  in  a  notional  form),  an  empirical  adjustment  technique  has  been  developed  to 
improve  noctiunal  temperature  gradient  estimates.  In  this  technique  the  parameter  is 
used,  where 

Qq  =  V^9P/P  -  ^zp/p.  {B  -  31) 

(Qq  is  the  right-hand  side  of  equation  (B-10).  This  parameter  is  then  used  in  the  equation 

X  =  min(1.5,  Qlidufdz)).  {B  -  32) 
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The  temperature  gradient  is  then  adjusted  (for  non-negative  temperature  gradients)  by 
the  equation 

^  =  ^(1 -0.125 X -1-0.376 X*).  (B-33) 

For  X  <  1/3,  the  maximum  difference  between  this  correction  and  the  normal  temperature 
gradient  calculation  occurs  at  X  =  1/6  where  the  correction  factor  is  1  —  1/96.  B^ond 
X  =  1/3  the  correction  factor  increases  to  the  mmfimnm  of  1.5  at  X  =  4/3.  This  value 
for  X  represents  a  point  1/3  beyond  the  X  :=  1  value  where  the  similarity  equations  fail. 
This  equation  therefore  induces  small  vmderestimating  error  for  X  <  1/3,  but  for  large 
X  the  error  between  measured  and  predicted  temperature  gradient  is  reduced.  Since 
the  correction  term  is  dependent  on  the  variable  that  influences  the  ability  of  turbulent 
turnover  of  small  turbules  as  discussed  in  section  5.  of  this  appendix,  the  correction  thus 
increases  the  temperature  gradient  to  compensate  for  the  error  induced  by  setting  u«  to  a 
minimum  value  to  avoid  negative  L  values. 
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